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ABSTRACT 


This  work  has  been  an  investigation  of  the  applica¬ 
bility  of  the  finite  element  method  to  the  study  of  crystal 
growth  processes.  The  finite  element  method  was  chosen 
because  a  review  of  past  analytical  research  in  the  area  of 
crystal  growth  theory  indicated  that  analytical  techniques 
may  have  reached  the  limit  of  their  usefulness  and  further 
progress  in  the  area  may  require  the  exploitation  of  numer¬ 
ical  methods. 

To  demonstrate  its  applicability,  the  method  has  been 
used  to  analyze  the  problems  such  as,  the  stability  of 
cylindrical  and  spherical  particles  growing  in  supercooled 
liquid,  dendritic  growth,  and  the  evolution  of  dendritic 
shape  with  time.  Quantitative  comparisons  between  this 
approach,  other  theories  and  experimental  data  are  made  for 
all  examples  to  assess  the  value  of  the  method. 

The  method  has  also  been  used  to  obtain  a  time  in¬ 
variant  shape  of  a  dendrite  growing  in  a  supercooled  liquid. 
The  results  indicate  that  the  velocity  of  growth  of  a  dendrite 
with  time  invariant  shape  is  closer  to  that  predicted  by  the 
Langer-Krumbhaar  stability  criterion  than  to  that  predicted  by 
the  maximum  velocity  principle  used  by  many  researchers. 

These  results  on  the  applicability  of  the  finite 
element  method  to  crystal  growth  problems  are  encouraging. 
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However,  certain  concerns  still  remain  regarding  accuracy 
in  calculations  of  surface  curvatures  and  separation  of 
numerical  instabilities  from  physical  instabilities.  If 
these  concerns  and  problems  can  be  overcome,  a  number  of 
interesting  applications  of  the  finite  element  method  to  the 
study  of  crystal  growth  phenomena  are  possible. 
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CHAPTER  1 


INTRODUCTION 

The  problems  related  to  the  crystallization  mechanism 
of  supercooled  liquids  have  been  studied  theoretically  and 
experimentally  for  some  time  by  physical  chemists  and 
materials  scientists.  In  recent  years,  however,  the  subject 
of  crystal  nucleation  and  growth  has  received  a  large  amount 
of  attention  from  scientists  and  engineers  working  in  such 
diverse  areas  as  solid  state  physics,  medicine,  biological 
sciences,  atmospheric  sciences,  solar  energy  utilization 
(thermal  energy  storage) ,  and  desalination  of  ocean  water. 

The  overall  crystallization  process  is  determined  by  transport 
effects  (heat  and  mass  transfer)  as  well  as  interface  effects 
(surface  tension  and  growth  kinetics) .  Discussions  of  mechan¬ 
ism  of  crystal  growth  must  consider  both,  as  well  as  the 
nucleation  process  itself.  As  Mott  [1]  pointed  out,  the 
transport  process  may  usually  be  treated  as  macroscopic, 
involving  the  processes  of  diffusion  of  heat  or  matter.  On 
the  other  hand,  the  interface  process  involves  the  detailed 
atomic  structure  of  the  surface. 

Depending  on  the  system  studied,  either  the  transport 
or  the  interface  effects  may  be  the  slower  and  hence  will  be  the  rate 
controlling  process*  Alternatively,  both  processes  may  control 
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the  rate  together  and  indeed,  interact  with  each  other  [2]. 

In  many  situations,  however,  the  rate  of  removal  of  latent 
heat  controls  the  rate  at  which  crystallization  can  proceed. 

If  this  is  the  case  the  problem  of  determining  the  rate  at 
which  the  interface  boundary  advances  through  the  medium 
reduces  to  one  of  calculating  the  temperature  distribution 
in  the  solid  and  liquid  phases. 

An  important  factor  controlling  the  micro-  and  macro¬ 
structure  of  the  crystal  grown  from  the  melt  is  the  direction 
of  the  heat  flux  at  the  interface  boundary.  For  solidification 
of  a  melt  which  is  above  its  equilibrium  fusion  temperature, 
heat  flows  from  the  melt  to  the  interface  and  solidification 
must  be  maintained  by  conduction  of  heat  from  the  interface 
into  the  crystal.  For  a  one  component  system  in  which  heat 
transfer  is  the  limiting  process  the  resulting  interface  tends 
to  be  smooth  and  isothermal.  However,  for  a  crystal  growing 
into  a  supercooled  melt  (temperature  less  than  the  equilibrium 
fusion  temperature)  heat  flows  from  the  interface  to  the  melt. 
In  this  case  the  growing  interface  will  tend  to  be  unstable 
and  dendritic  growth  will  result.  It  is  to  problems  of  in¬ 
stability  and  dendritic  growth  that  this  thesis  will  be 
addressed . 

The  term  'dendritic  growth'  when  used  in  reference  to 
solidification  suggests  a  tree-shaped  crystal  structure,  that 
is  a  crystal  with  an  intricate  network  of  branches.  Dendritic 
solidification  is  characterized  by  a  morphology  resulting 
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from  the  growth  of  long,  thin  spikes  in  specific  crystallo¬ 
graphic  directions,  with  regular  branches  in  other  equivalent 
directions.  The  branching  habit  extends  to  secondary,  tertiary, 
and  sometimes  higher  orders.  The  main  qualitative  experiment¬ 
al  observations  are:  (a)  that  for  a  pure  melt,  dendritic  growth 
takes  place  only  when  the  melt  is  supercooled,  (b)  that  the 
directions  of  growth  are  always  strictly  crystallographic,  (c) 
branching  occurs  at  roughly  regular  spacing,  smaller  for  each 
successive  order  of  branching,  and  (d)  that  only  a  small  portion 
of  liquid  solidifies  in  this  way  [3] .  Crystal  growth  from 
all  three  states  of  matter — vapour ,  liquid  and  solid — is  known 
to  produce  this  type  of  morphology. 

The  subject  of  dendritic  growth  has  a  whole  series  of 
important  consequences  for  the  growth  of  crystals  of  many 
commonly  used  materials.  Some  idea  of  the  scale  of  this  can 
be  judged  from  the  fact  that  hundreds  of  million  tons  of 
steel  are  produced  annually,  all  of  which  has  undergone 
dendritic  growth.  There  is  a  wide  range  of  data  available 
in  the  literature  on  the  growth  of  dendrites  in  different 
systems  (metals,  ice,  phosphorus,  solid  state  precipitates, 
and  ionic  crystals  from  aqueous  solutions) .  These  data  give 
an  understanding,  in  a  qualitative  sense,  of  the  basic  pheno¬ 
mena  associated  with  dendritic  growth.  However,  to  examine 
this  understanding  from  a  quantitative  viewpoint,  the  complex 
heat  flow  problem  with  moving  boundary  has  to  be  solved.  In 
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a  general  theoretical  examination  of  the  problem  the  essential 
difficulty  resides  in  the  fact  that  the  position  of  the 
boundary  on  which  one  wishes  to  apply  boundary  conditions 
is  not  known;  finding  its  position  is  part  of  the  problem. 

This  makes  the  problem  non-linear  so  that  solutions  may  not 
be  built  up  by  superposition.  Consequently,  analytical 
solutions  must  be  sought  through  application  of  approximate 
methods . 

A  review  of  past  analytical  research  in  the  area  of 
dendritic  growth  theory  suggests  that  these  techniques  may 
have  reached  the  limit  of  their  usefulness  and  further  pro¬ 
gress  in  the  area  may  require  the  exploitation  of  numerical 
methods.  It  is  believed  that  the  Finite  Element  Method  is 
well  suited  to  solve  such  problems,  since  the  basic  concepts 
of  this  method  have  already  been  found  to  possess  general 
applicability  to  a  wide  range  of  field  problems.  Indeed, 
in  addition  to  the  usual  applications  in  structural  and 
continuum  mechanics  [4],  use  can  also  be  found  in  many  other 
fields.  Problems  in  such  diverse  fields  as  water  seepage  [5], 
irrotational  flow  of  ideal  fluids  [6],  electrical  fields  [7], 
and  heat  conduction  [8]  have  already  been  studied  by  the 
finite  element  approach. 

It  is  the  purpose  of  the  present  work  to  examine  the 
applicability  of  the  finite  element  method  to  the  study  of 
crystal  growth  transport  processes.  Specifically,  the  method 
will  be  tested  on  two  types  of  problems.  The  first  is  an 
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investigation  of  the  morphological  stability  of  cylindrical 
and  spherical  particles  undergoing  radial  growth  controlled 
by  diffusion.  The  second  is  the  prediction  of  the  growth 
velocity  of  a  fully  developed  dendrite  as  a  function  of 
supercooling.  The  primary  motivation  is  to  see  if  the  finite 
element  method  can  improve  our  understanding  of  the  crystalli¬ 
zation  process. 

To  this  end,  a  brief  literature  survey  is  made  in 
Chapter  2  in  order  to  establish  a  relation  between  the  present 
work  and  what  has  been  reported  previously.  Chapter  3  is 
devoted  to  fundamental  theory,  which  includes  definitions 
and  basic  equations  relevant  to  this  study,  the  variational 
principle,  and  fundamental  concepts  of  finite  element  analysis 
and  the  Ritz  technique.  Chapter  4  presents  the  analysis  of 
morphological  stability  of  spherical  and  cylindrical  particles 
growing  in  supercooled  liquids.  Chapter  5  contains  the  analy¬ 
sis  of  dendrites  with  the  shapes  of  parabolic  cylinders  and 
paraboloids  of  revolution.  In  Chapter  6,  an  attempt  has  been 
made  to  demonstrate  that  the  finite  element  method  can  be  used 
to  study  the  time  dependent  growth  of  dendrites. 

The  results  obtained  by  the  present  approach  are 
compared  with  those  previously  obtained  by  other  methods, 
either  theoretical  or  experimental.  It  is  found  that  good 
agreement  exists  among  the  results  predicted  by  these  different 
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Based  on  the  fact  that  the  results  compare  well 
for  the  problems  that  were  treated,  it  is  suggested  that 
the  finite  element  method  may  be  developed  into  a  useful 
tool  in  the  analysis  of  other  crystal  growth  problems. 


CHAPTER  2 


LITERATURE  SURVEY 

2 . 1  Introduction 

A  review  of  the  most  significant  contributions  to  the 
area  of  morphological  stability  in  crystal  growth  and  in 
particular  to  that  of  dendritic  growth  will  be  presented  here. 
Due  to  the  vast  literature  associated  with  crystal  growth 
problems,  only  those  papers  which  are  pertinent  to  the 
present  investigation  will  be  cited.  An  excellent  survey  of 
literature  on  crystal  growth  mechanisms  is  given  in  the  works 
of  Parker  [2]  and  Jackson  [9].  For  a  historical  review  of 
the  general  problem  of  dendritic  growth,  the  reader  is  referred 
to  the  study  by  Smith  [10] . 

2 . 2  Experimental  Investigations 

The  term  dendrite  was  apparently  first  introduced  to 
the  world  of  crystal  growth  by  Tschernoff  in  1879  [10].  He 
used  it  to  describe  the  highly  branched  structure  he  found 
in  the  centre  of  metal  ingots.  With  the  recognition  of  crystal 
structure  of  metals  it  was  assumed,  and  later  demonstrated, 
that  the  directions  of  the  various  dendrite  arms  coincided 
with  crystographically  significant  directions  of  the  structure. 
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The  modern  impetus  to  this  area  of  crystal  growth 
science  may  be  traced  to  the  experimental  studies  of  Chalmers 
and  his  collaborators  [11-16] ,  who  demonstrated  that  the 
dendrites  represent  the  most  advanced  stage  of  interfacial 
instabilities  in  a  wide  class  of  materials.  With  the  help  of 
experiments,  Wineburg  and  Chalmers  [11,12]  also  showed  that 
dendritic  growth  is  observed  only  if  the  liquid  is  at  a  lower 
temperature  than  the  interface. 

2.2.1  Dendritic  Growth  Rate 

The  growth  rate  of  dendrites  for  various  amounts  of 
supercoolings  has  been  measured  in  pure  water  [17],  tin  [14,18], 
nickel  [19],  lead  [14],  and  plastic  crystals  [20].  However, 
some  of  these  investigations  dealt  with  the  growth  rate  of 
crystals  in  contact  with  a  solid  substrate,  which  bear  little 
relationship  to  the  results  obtained  in  free  growth  from  the 
melt . 

Wineburg  and  Chalmers  [12]  measured  the  free  growth 
of  lead  dendrites,  and  showed  that  the  dendrites  grew  much 
faster  than  the  smooth  interface.  Rosenburg  and  Winegard 
[18]  measured  the  growth  rate  of  dendrites  in  a  supercooled 
bath  of  tin  with  a  range  of  supercooling  from  0.4°C  to  11°C. 
They  found  that  there  was  a  substantial  scatter  in  the  results, 
and  a  factor  of  three  existed  between  the  highest  and  lowest 
values  for  a  given  temperature. 
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Lindenmeyer  and  Chalmers  [13]  measured  the  free  growth 
of  ice  in  pure  water  along  the  basal  plane  at  bath  super¬ 
coolings  between  2°'  and  6.5°C.  They  reported  the  free  growth 
rate  v  in  distilled  water  to  be 

v  =  0.028  (AT)2*39  cm/s ec,  (2-1) 

where  AT  is  in  °C.  In  distilled  water  having  similar  super¬ 
cooling,  they  observed  that  the  growth  rate  upon  a  substrate 
was  as  much  as  an  order  of  magnitude  faster.  It  was  surmized 
that  this  difference  was  due  to  the  proximity  of  a  heat  sink 
which  allowed  a  more  efficient  removal  of  latent  heat.  They 
concluded  from  this  result  that  the  rate  of  growth  of  ice 
crystals  in  water  is  limited  by  the  conduction  of  latent  heat 
into  the  water  for  free  growth,  and  into  the  water,  the  solid 
wall,  and  across  the  solid  wall  into  the  bath  for  growth 
on  a  substrate. 

Walker  [19]  made  extensive  measurements  on  the  growth 

Of  dendrites  in  nickel  and  cobalt,  and  plotted  the  growth 

velocity  against  the  square  of  the  supercooling.  He  found  that 

up  to  a  supercooling  of  almost  175 °C,  the  results  fell  on 

2 

the  v  «  (AT)  line. 

Ryan  [21]  investigated  the  free  growth  of  single  ice 
crystals  in  supercooled  water,  and  measured  growth  rates  for 
supercoolings  ranging  from  0.1  to  2.5°C.  His  data  do  not 
conform  to  the  v  «  (AT)n  law  for  the  entire  range  of  super¬ 
cooling.  In  the  range  of  0.1°C  to  1.3°C,  the  growth  rate  obeys 
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approximately  a  power  law  of  supercooling.  Ryan  also  reported 
that  in  the  range  of  subcoolings  from  0.2  to  1.5°C  the  shape 
of  the  dendrite  tip  may  be  approximated  by  a  parabolic 
platelet.  His  results  also  showed  that  the  tip  radius  de¬ 
creases  as  the  supercooling  increases. 

Mason  [22]  has  reported  that  ice  crystals  grow  as 
needles  (paraboloids  of  revolution)  when  grown  in  water  in 
the  supercooling  range  from  3  to  5°C.  Hallet  [23]  measured 
the  growth  rates  of  ice  dendrites  growing  parallel  to  the 
basal  plane  in  supercooled  water,  and  his  data  followed  the 
relation 

v  =  0.08  (AT)1*9  cm/sec,  (2-2) 

for  supercooling  between  0.1° C  and  2.0°C.  He  also  found 
that  the  dendrite  arm  width  and  spacing  decreased  with  in¬ 
creased  supercoolings. 

Macklin  and  Ryan  [24]  made  a  detailed  study  of  the 
growth  of  ice  in  pure  water  at  supercoolings  up  to  7.5°C. 

They  reported  growth  velocities  which  approximately  fit  the 
relationship 

v  =  0.0227  (AT)2*29  cm/sec ,  (2-3) 

where  AT  is  in  0.°C. 

Pruppacher  [25]  determined  the  rate  at  which  ice 
crystals  grow  freely  in  supercooled  water  and  in  dilute 
acqueous  solutions  of  various  salts.  He  reported  that  at 
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supercoolings  between  0.5°C  and  9°C,  the  growth  rate  of  ice 
dendrites  can  be  represented  by  the  relation 

v  =  0.035  (AT)2*22  cm/sec,  (2-4) 

He  also  concluded  that  dissolved  salts  affected  the  growth 

rate  of  ice  crystals  in  a  systematic  manner.  At  concentrations 

-2 

larger  than  10  moles/liter  all  salts  reduced  the  growth 

rate  below  the  one  in  pure  water  by  an  amount  which  depended 

on  type  of  salt  in  solution.  At  concentrations  smaller  than 
-2  .  . 

10  moles/liter  some  salts  did  not  affect  the  growth  rate 
whereas  others  increased  it  over  and  above  in  the  pure  water. 

He  concluded  that  at  low  solute  concentrations  the  growth 
rate  of  ice  in  aqueous  solutions  is  limited  by  the  rate  of 
dissipation  of  latent  heat.  As  the  solute  concentration 
increased,  the  rate  of  diffusion  of  solute  atoms  away  from  the 
interface  became  increasingly  important. 

Glicksman  and  Schaefer  [26]  have  reported  some 
results  on  the  in  situ  observations  of  dendritic  solidification 
in  pure  tin.  They  used  cinematographic  techniques  to  obtain 
quantitative  data  on  the  variation  of  tip  profile  as  a  function 
of  the  growth  speed,  and,  on  the  transition  of  periodic  edge- 
profile  perturbations  into  true  dendritic  side  branches.  They 
observed  that  the  approximation  of  a  parabolic  dendrite  tip 
appeared  satisfactory  over  a  distance  of  about  9  tip-radii 
back  from  the  nose  of  each  dendrite.  Further  back  than  about 
9  tip-radii  the  dendrite's  shape  became  increasingly  distorted 
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from  a  parabola  by  the  combined  influences  of  side-branching 
and  interferring  neighboring  thermal  diffusion  fields. 

Kallungal  [27]  and  Kallungal  and  Barduhn  [28]  con¬ 
ducted  extensive  experiments  to  measure  the  growth  rate  of 
ice  crystals  in  slightly  supercooled  water.  They  observed 
that  the  shape  of  ice  dendrites  was  neither  2-D  parabolic 
cylinder  nor  axisymmetric  paraboloid  of  revolution,  but  it 
was  similar  to  an  elliptic  paraboloid.  They  also  reported 
that  the  ice  crystal  growth  in  slightly  subcooled  quiscent 
water  could  be  represented  by  the  relation 

v  =  0.0118  (AT)2*17  cm/sec,  (2-5) 

where  AT  is  in  °C. 

Kallungal  and  Barduhn  [28]  also  measured  the  tip 
radius  of  curvature  of  ice  crystals  and  observed  that  the 
tip  radii  could  be  expressed  by  the  relation 

R  =  0.61/AT  microns  (2-6) 

where  AT  is  in  °C  and  is  in  the  range  of  0.1°  to  1°C. 

Glicksman,  Schaefer,  and  Ayers  [20]  designed  an 
experimental  set  up  to  study  the  dendritic  growth  in  super¬ 
cooled  succinonitrile ,  a  low  entropy-of-fusion  plastic  crystal. 
They  measured  the  growth  rate  of  dendrites  for  free  growth 
conditions  as  well  as  confined  growth  conditions.  They 
reported  that  for  free  growth,  the  growth  rate  data  fitted 


a  relation 


V  a  (AT) 


2.6 


cm/sec , 


(2-7) 


for  supercoolings  ranging  from  0.1°C  to  7°C.  They  also 
measured  the  tip  radii  of  curvature,  and  reported  that  the 
tip  radius  of  curvature  was  inversely  proportional  to  the 
amount  of  supercooling.  Furthermore,  they  found  that  near 
the  tip,  the  dendrites  were  almost  paraboloids  of  revolution. 
They  compared  their  experimental  results  with  three  different 
steady-state  theories  and  found  that  full  agreement  between 
steady-state  theory  and  experiment  was  clearly  lacking.  They 
reported  that  the  velocities  obtained  by  theoretical  methods 
were  considerably  higher  than  the  experimental  values.  How¬ 
ever,  they  found  that  all  the  non-isothermal  steady-state 
theories  considered,  predicted  the  correct  Peclet  number¬ 
supercooling  relationship  over  a  limited  range  of  supercooling. 
They  also  reported  that  the  interfacial  molecular  attachment 
kinetics  are  sufficiently  rapid  not  to  influence  significantly 
the  axial  growth  rate  of  freely  growing  dendrites  of  succino- 
nitrile  over  a  limited  range  of  supercooling. 

Summarizing  the  experimental  work  on  the  dendritic 
growth  rate  it  can  be  concluded  that  although  no  two  investi¬ 
gators  agree  on  the  reported  value  of  the  growth  rates ,  the 
experimental  data  usually  conform  to  a  curve  of  the  form 

v  =  A  ATn  (2-8) 

where  v  is  the  growth  velocity,  AT  is  the  supercooling  and 
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A  and  n  are  constants  for  a  given  material.  For  all  the 
measurements,  n  lies  between  1.5  and  3.0.  Also,  the  shape  of  the 
dendrite  tip  is  somewhere  between  a  paraboloid  of  revolution 
and  a  2-D  parabolic  cylinder  depending  on  the  material.  Finally, 
the  radius  of  curvature  at  the  tip  of  dendrites  is  reported 
to  be  proportional  to  1/AT  by  all  the  investigators. 

2.2.2  Morphological  Stability 

Morphological  stability  in  crystal  growth  is  concerned 
with  the  stability  of  a  particular  shape  of  solid,  growing 
from  melt,  solution,  or  vapour  by  diffusional  process,  that 
is,  whether  it  is  stable  against  small  perturbations  in  shape. 

Arakawa  [29],  and  Arakawa  and  Higuchi  [30]  studied 
the  growth  of  ice  crystals  in  water  and  observed  that  when  the 
supercooled  water  (temperature  about  -0.5°C)  was  nucleated 
by  placing  small  hoar  frost  crystals  on  the  surface,  needle 
and  circular  disc  crystals  grew  outwards  from  the  nucleating 
crystals.  The  circular  disc  crystals  developed  notches 
around  their  edges  as  they  grew,  and  later  they  became  six- 
sided  dendrites. 

Williamson  and  Chalmers  [31]  photographed  growing 
ice  discs  in  slightly  subcooled  water  and  observed  that  at 
subcoolings  greater  than  0.4°C,  the  disc  morphology  became 
unstable  and  small  protuberances  appeared  on  the  edge  of  the 
disc.  At  subcoolings  greater  than  0.6°C/ice  assumed  the 
characteristic  dendritic  appearance  shown  by  snow  flakes. 
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They  reported  that  as  the  supercooling  was  made  greater 
than  1°C,  dendrites  with  side-branches  appeared  in  specifi¬ 
cally  preferred  growth  directions.  Also,  as  the  subcooling 
was  increased  the  spacing  between  dendrite  arms  decreased 
and  they  found  that  the  dendrite  arm  width  and  spacing 
obeyed  the  following  relation  approximately: 

-2 

arm  width  10  /AT  cm 

-2 

arm  spacing  10  /AT  cm 

Lindenmeyer  and  Chalmers  [32]  showed  that  a  circular 
crystal  became  hexagonal  in  shape  and  finally  developed  into 
a  dendrite  crystal  with  six  arms. 

Morris  and  Winegard  [33]  photographed  the  growing 
dendrites  and  observed  that  the  dendrites  tip  grew  with  a 
variable  shape,  and  fluctuated  in  a  periodic  manner  as  it 
propagated  into  the  melt.  From  their  photographs  it  is 
evident  that  the  tip  continuously  grew  larger  in  radius  and 
subsequently  small  irregularities  occurred  which  eventually 
grew  into  primary  branches. 

Hardy  and  Coriell  [34,35]  have  made  measurements  of 
the  growth  (and  decay)  of  perturbations  on  single  crystal 
clinders  of  ice  growing  from  slightly  supercooled  distilled 
water.  Ice  crystals  were  grown  with  the  basal  plane  normal 
to  the  axis  of  the  cylinder.  Hence,  growth  was  occurring 
parallel  to  this  plane;  the  kinetic  coefficient  for  such 


16 


growth  is  known  to  be  sufficiently  large  so  that  the  kinetic 
(interface  attachment)  contribution  to  the  undercooling  was 
negligible  in  those  experiments.  Due  to  the  low  undercoolings 
and  small  growth  velocities,  the  conditions  were  close  to 
steady-state.  The  ice  crystal  was  first  allowed  to  come  into 
equilibrium  with  a  large  bath  of  water  at  a  temperature 
slightly  above  0°C.  The  temperature  of  the  bath  was  decreased 
by  about  0.1°C  and  the  ice  crystal  started  to  grow  from  radius 
of  about  0.05  cm.  The  crystal  could  be  rotated  about  its 
axis,  and  photographed  from  the  side  so  the  radius  and  size 
of  the  perturbations  could  be  measured.  They  reported  that 
at  first  the  ice  crystal  grew  as  a  regular  cylinder  and  the 
solid-liquid  interface  could  be  characterized  by  an  equation 
(in  cylindrical  coordinates  r,(j),z)  of  the  form  r  =  R  where  R 
depends  on  time  but  not  on  z  and  <J> .  As  the  cylinder  grew, 
instabilities  developed  and  the  topography  of  the  unstable 
interface  could  be  approximately  characterized  by  the  equation 

r  =  R  +  6  cos(KcJ))  (2-9) 

where  K=6  or  12,  and  6  is  a  small  amplitude  factor.  At  a 
later  stage, perturbations  also  developed  along  the  z  direction. 

In  concluding  the  experimental  literature  survey,  it 
should  be  pointed  out  that  though  a  wide  range  of  experimental 
data  exist  on  the  problems  of  morphological  stability  and 
dendritic  growth,  care  must  be  exercised  when  comparing  them 
with  theory  because  the  precise  conditions  under  which  the 


experiments  were  conducted  are  often  not  known. 
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2 . 3  Theoretical  Investigations 

While  the  bulk  of  experimental  investigations  of 
dendritic  growth  have  been  about  equally  divided  between 
free  growth  and  confined  growth,  this  certainly  is  not  true 
of  theoretical  investigations  reported  in  literature  since 
the  overwhelming  majority  of  these  works  have  dealt  with 
free  growth  conditions.  As  mentioned  earlier,  either  diffu¬ 
sion  process  or  interface  kinetics  limit  the  growth  rate 
under  free  growth  conditions.  For  this  paper  we  will  review 
only  those  theories  which  consider  the  diffusion  process  as 
the  rate  controlling  process  as  the  crystallization  proceeds, 
although  we  will  mention  the  use  of  interface  kinetic  laws  as 
boundary  conditions  for  the  transport  process. 

2.3.1  Dendritic  Growth  Rate 

Papapetrou  [36]/  in  a  classical  study  of  dendritic 
growth,  suggested  that  the  tip  of  a  dendrite  has  the  form  of 
a  paraboloid  of  revolution.  But  it  was  G.P.  Ivantsov  [37] 
of  the  Soviet  Union  who  made  the  first  serious  attempt  to  make 
a  mathematical  analysis  of  dendritic  growth  rates  into  its 
supercooled  melt.  He  showed  that  some  aspects  of  dendritic 
growth  could  be  described  quantitatively  with  heat  flow  theory 
applied  at  a  suitable  microscopic  level.  Using  the  suggestion 
of  Papapetrou  [36],  he  assumed  that  the  tip  region  of  a 
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dendritic  trunk  may  be  adequately  represented  by  a  paraboloid 
of  revolution.  He  further  assumed  that  the  crystal  was 
isothermal  and  advancing  at  a  constant  velocity  v,  into  an 
isothermal,  undercooled  bath  at  temperature  T  .  He  obtained 
the  solution  to  the  moving  boundary  problem  diffusion  equation 
and  showed  that  the  dendrite  obeys  the  equation 

-vR/2a£  exp  (vR/2a £)  Ei(-vR/2a£)  =  AT  c£/L  (2-10) 

where  v  is  the  growth  rate,  R  is  the  dendrite  tip  radius, 
a^  and  c^  are  the  thermal  diffusivity  and  the  specific  heat 
of  the  melt,  L  is  the  latent  heat  of  fusion,  AT  is  the  amount 
of  bath  supercooling,  and  is  the  integral  error  function. 

However,  the  above  solution  does  not  completely  specify  the 
conditions  of  growth.  A  diffusion  solution  can  be  obtained 
for  a  given  undercooling  provided  the  product  vR  is  constant, 
as  is  evident  from  above  equation  (2-10).  Thus,  for  a  given 
bath  supercooling,  thin  dendrites  growing  rapidly  or  thick 
dendrites  growing  slowly  can  both  satisfy  the  diffusion 
conditions.  Commonly,  a  value  of  R  which  maximizes  v  would 
be  sought  but  this  is  impossible  since  vR  =  constant. 

In  1960,  D.  E.  Temkin  [38],  also  of  Soviet  Union, 
recognized  the  limitations  of  the  Ivantsov'  analysis  and  argued 
that  the  assumption  of  an  isothermal  interface  was  fundamentally 
incorrect.  He  noted  that  the  very  existence  of  the  variation 
in  curvature  along  the  surface  of  the  dendrite  implies  that 
the  equilibrium  fusion  temperature  must  vary  with  position  in 
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accordance  with  the  Gibbs-Thompson  relationship.  Temkin 
also  introduced  linear  interface  kinetics  into  the  problem. 
For  mathematical  simplicity,  Temkin  assumed  that  the  shape 
is  still  a  paraboloid  of  revolution.  He  also  assumed  that 
Laplace  equation  could  be  used  for  the  problem  with  small 
bath  undercoolings,  that  is  when  vR/a^  <<  1.  He  obtained 
the  solution  to  the  diffusion  problem  with  above  mentioned 
condition  to  be 


vR/2a^  exp  (vR/2a^)  E^(-vR/2a^) 


AT 

L  p 

s 


_ 1 

l+L1/yR+L2y/vR2 


(2-11) 


where  y  is  the  interface  kinetic  coefficient,  y  is  the 
specific  surface  free  energy  and  and  pg  are  the  density 

of  the  liquid  and  the  solid  phase  respectively.  The  coeffi¬ 
cients  and  L2  are  defined  by 
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where  T^  is  the  fusion  temperature  of  the  pure  bulk  material 

and  k  and  k  are  the  thermal  conductivities  of  the  liquid 
l  s 


and  solid. 
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The  major  difference  between  the  Temkin  and  Ivantsov 
results  is  that  the  former  exhibits  a  maximum  in  the  function 
v  versus  R.  Temkin  chose  the  value  of  R  to  be  the  actual 
tip  radius  which  corresponded  to  the  maximum  velocity  of 
growth.  For  vR/2a^  <<  1,  Temkin  showed  that  the  equation 
(2-11)  may  be  written  in  the  following  approximate  form, 


vR 

2a 


a 


^  AT  VL _ 

Ps  (1  +  I^/yR  +  L2  y/vR2) 


(2-14) 


where  =  0.457,  and  n  =  1.21.  Temkin  applied  the  maximum 
growth  rate  principle  to  equation  (2-14)  to  get, 


( 2n-l )  L2  y 

v  =  - - - 

R  [R  -  (n-1)  L1/y] 


(2-15) 


and  R  is  obtained  from  the  relation 
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Thus,  according  to  Temkin,  the  simultaneous  solution 
of  equations  (2-15)  and  (2-16)  uniquely  predicts  both  the 
velocity  and  tip  radius  of  a  freely  growing  dendrite  for  a 
given  bath  undercooling . 

For  the  special  case  where  the  kinetic  coefficient, 
y,  is  very  large  we  may  neglect  the  term  L^/y  in  equations 


-**  J 
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(2-15)  and  (2-16)  to  obtain 


(2n-l)  L2y/v, 
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(2-18) 


At  approximately  the  same  time  as  the  Temkin  paper, 
Bolling  and  Tiller  [16]  presented  an  excellent  physical 
picture  of  the  effects  of  non-isothermality  on  the  growth 
of  a  dendrite.  Independent  of  Temkin  they  also  recognized 
the  non-isothermality  of  a  dendrite  due  to  the  variation  of 
curvature  and  kinetics  over  the  surface.  They  showed  that 
the  variation  of  temperature  in  the  solid  gives  rise  to  a 
heat  flux  in  addition  to  the  latent  heat  term  at  the  tip  of 
the  dendrite.  Thus,  the  velocity  of  a  non-isothermal  dendrite 
will  be  reduced  compared  to  an  isothermal  dendrite  under 
similar  conditions.  Bolling  and  Tiller  also  attempted  to 
derive  expressions  which  would  determine  both  the  maximum 
velocity  and  the  corresponding  tip  radius.  Inspite  of  the 
fact  that  in  a  later  analysis  Kotler  and  Tarshis  [39]  have 
shown  that  their  calculations  had  some  mathematical  inconsis¬ 
tencies,  the  Bolling  and  Tiller  [16]  work  contains  numerous 
useful  contributions  to  the  understanding  of  dendritic  growth. 
Bolling  and  Tiller  gave  a  detailed  review  and  discussion  of 
the  following  general  assumptions  made  in  the  analysis  of 
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dendritic  growth  by  theoretical  investigators  in  this  field. 
These  assumptions  are: 

(a)  growth  occurs  at  a  constant  speed  (quasi¬ 
stationary  process) , 

(b)  the  dendrite  shape  is  either  a  parabolic  cylinder 
or  paraboloid  of  revolution, 

(c)  the  liquid  domain  is  infinite, 

(d)  the  maximum  velocity  condition  is  used  to 
specify  the  velocity, 

(e)  natural  convection  due  to  density  differences 
in  solid  and  liquid  phase  is  neglected. 

Horvay  and  Cahn  [40]  analyzed  the  dendrite  growth 
problem  using  four  different  geometries  for  the  crystal  shape; 
namely  (a)  parabolic  cylinder,  (b)  elliptic  paraboloid,  (c) 
paraboloid  of  revolution,  (d)  spheroidal.  However,  they  consid¬ 
ered  the  dendritic  interface  to  be  isothermal. 

Treating  the  problem  in  parabolic  coordinates  in  a  moving 
frame  of  reference;  they  used  the  conventional  technique  of 
separation  of  variables.  For  each  geometry,  they  were  able 
to  express  the  growth  rate  as 

v  =  a  (AT)  n/R  cm/sec .  (2-19) 

o 

where  AT  is  in  °C.  The  factor  a  varied  by  almost  two  orders 

o 

of  magnitude  between  the  parabolic  cylinder  and  the  paraboloid 


of  revolution. 
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Hillig  [41]  solved  the  steady-state  growth  of  a 


dendrite  platelet  into  its  pure  supercooled  melt  by  using 
what  he  called  a  "self-consistent"  method.  This  method 
requires  solution  for  the  shape  of  the  crystal  that  simultan¬ 
eously  satisfies  the  required  heat  flow  and  temperature 
variation  over  the  dendrite  surface.  In  his  method,  the 
temperature  at  infinity  was  not  assumed  but  calculated 
using  a  Taylor  expansion.  He  concluded  that  the  assumption 
of  steady  state  is  probably  invalid.  He  suggested  that 
possibly  the  shape  is  oscillatory  in  time  and  the  dendrite 
moves  with  a  constant  average  velocity. 


Trivedi  [42]  presented  an  exact  solution  to  the  heat 


flow  equation  for  the  growth  of  dendritic  needles  (paraboloids 
of  revolution)  in  a  supercooled  melt.  His  final  result  is 
similar  to  the  equation  [2-11],  but  his  values  for  and  L 2 
are  not  constants  but  functions  of  the  thermal  Peclet  number. 
He  also  compared  his  theoretical  results  with  experimental 
results  on  phosphorus  and  obtained  good  agreement. 


Trivedi  [43]  also  analyzed  the  growth  of  non-isothermal 


dendritic  plates  (parabolic  cylinders)  growing  in  a  super¬ 
cooled  melt.  He  expressed  the  solution  for  this  problem  as 


1 


L 


(2-20) 


where  1^  and  L2  are  of  the  form  defined  by  equation  (2-12) 

and  equation  (2-13) ,  but  have  different  values  from  the  case 

of  the  paraboloid  of  revolution.  For  the  case  where  kinetic 

effects  may  be  neglected  and  the  Peclet  number  is  small, 

2a£ 

the  above  result  can  be  simplified  to  give 
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Differentiating  equation  (2-21)  with  respect  to  R  and  setting 
dv 

dR  =  0  ^1VeS 

2 

R  =  3yL2/v.  (2-22) 


Substituting  equation  (2-22)  into  equation  (2-21)  results  in 
the  growth  model  for  a  parabolic  cylinder  as 
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Holzman  [44,45,46]  has  made  analysis  similar  to  but 
independent  of  Trivedi  for  both  parabolic  platelet  and  para¬ 
boloid  of  revolution.  He  questioned  the  traditional  time 
invariant  shape  of  the  dendrites.  To  maintain  a  time  invariant 
shape,  the  heat  balance  at  the  solid-liquid  interface  requires 
that  the  sum  of  heat  fluxes  in  the  solid  and  liquid  at  each 
point  of  the  advancing  interface  must  be  equal  to  the  local 
rate  of  release  of  heat  of  fusion.  Holzman  found  that  this 
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is  not  true ,  and  he  calculated  a  quantity  he  called  "excess 
velocity" ,  which  is  the  difference  between  the  true  local 
growth  velocity  required  to  satisfy  heat  balance,  and  the 
growth  velocity  predicted  for  the  hypothetical  steady  state 
(shape  preserving)  model.  His  results  show  that  the  hypo¬ 
thetical  parabolic  or  paraboloidal  interface  tends  to  bulge 
with  a  pronounced  peak  at  about  one  radius  back  from  the 
tip  of  the  dendrite. 

The  analyses  of  dendritic  growth  which  we  have  just 
reviewed  have  all  been  mathematical  and  based  on  idealized, 
though  more  and  more  realistic  formulations  of  the  problem. 

With  the  advent  of  high  speed  computers  the  numerical  analysis 
of  more  realistic  models  became  possible.  Geering,  Oldfield 
and  Tiller  [47]  were  the  first  to  analyze  the  dendritic  growth 
problem  using  a  computer  model.  They  simplified  the  problem 
considerably  by  considering  only  one-dimensional  heat  transfer. 
This  was  done  by  choosing  the  subdivisions  (or  boxes)  of  the 
domain  so  that  two  sides  were  always  isotherms  and  the  other 
four  were  orthogonal  surfaces.  The  computer  algorithm  for 
choosing  the  boxes  in  this  way  was  analogous  to  seeking  an 
orthogonal  set  in  a  mathematical  solution.  They  checked  the 
results  obtained  from  the  computer  model  against  Ivantsov's 
solution  and  found  to  be  in  excellent  agreement.  They  also 
used  the  model  with  heat  flow  and  interface  curvature  effect 
taken  into  account  and  found  in  agreement  with  Temkin's 
analysis  for  pure  tin. 
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Oldfield  [48]  extended  the  above  mentioned  computer 
model  to  simulate  the  long  term  growth  of  tin  dendrites.  He 
obtained  the  smoothed  velocity  of  the  dendrite  tip  and  found 
that  the  growth  of  the  dendrite  was  accompanied  by  changes 
in  shape  and  a  fall  in  the  growth  velocity.  Oldfield  surmized 
that  the  dendrite  grows  with  a  cyclic  flucuation  in  velocity, 
forming  a  new  branch  in  each  phase  of  slow  growth.  However, 
the  above  mentioned  computer  models  were  quite  simplified  in 
form  and  the  heat  flux  in  the  solid  was  neglected. 

Recently,  Langer  and  Muller-Krumbhaar  [49-51]  have 

suggested  replacing  the  widely  used  maximum  velocity  principle 

.  .  .  .  2 
by  a  stability  criterion  of  the  form  vR  =  constant,  where 

v  is  the  growth  velocity  and  R  is  the  corresponding  tip  radius 

of  the  dendrite.  They  claim  that  this  type  of  stability 

criterion  is  valid  for  small  Peclet  numbers.  They  found  that 

a  direct  evaluation  of  the  relevant  constant  of  proportionality 

yields  values  of  v  and  R,  as  function  of  undercooling  which 

are  in  substantial  agreement  with  the  measurements  of  Glicksman 

et  al.  [20]  on  dendritic  growth  of  succinonitrile  crystals. 

Using  the  Langer-Muller-Krumbhaar  stability  theory  Langer, 

Sekerka  and  Fujioka  [52]  have  presented  a  universal  law  for 

dendritic  growth  rates.  They  checked  their  universal  law  with 

measured  values  on  ice  and  succinonitrile  dendrites  and  obtained 

good  agreement. 

Summarizing,  we  find  that  both  isothermal  and  non- 
isothermal  theories  predict  that  the  growth  rate  decreases  as 
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we  change  the  shape  from  a  paraboloid  of  revolution  to  that 


of  a  parabolic  cylinder.  The  exponents  on  subcoolings,  AT 
varies  from  2  to  4  and  the  growth  rates  predicted  at  a  given 
AT  vary  by  about  5  orders  of  magnitude. 

2.3.2  Morphological  Stability 

The  stability  is  tested  in  theoretical  development, 
by  assuming  an  arbitrarily  small  perturbation  in  shape  to  be 
present,  and  then  seeing  whether  this  perturbation  will  grow 
or  decay.  This  is  different  from  the  question  of  whether  the 
shape  and  equations  are  shape  preserving  (dendritic  growth 
problem) . 

Mullins  and  Sekerka  [53]  in  a  classical  paper,  intro¬ 
duced  into  crystal  growth  studies  the  concept  of  quantitative 
testing  for  shape  stability.  Prior  to  this  paper,  the  shape 
preserving  ellipsoids  of  Ham  [54]  were  generally  believed  to 
imply  that  such  shapes  were  stable  whereas  in  fact  they  had 
not  been  tested  for  stability. 

Mullins  and  Sekerka  [53]  studied  the  stability  of  the 
shape  of  a  spherical  particle  undergoing  diffusion  controlled 
growth  into  an  initially  uniformly  supersaturated  matrix. 

This  was  done  by  supposing  an  expansion,  into  spherical 
harmonics,  of  an  infinitesimal  deviation  of  the  particle  from 
sphericity  and  then  calculating  the  time  dependence  of  the 
coefficient  of  the  expansion.  It  was  assumed  that  the  pertin¬ 
ent  diffusion  field  obeyed  the  Laplace's  equation,  an  assumption 
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whose  conditions  of  validity  were  discussed  in  detail.  It 
was  shown  that  the  sphere  is  stable  below  and  unstable  above 
a  certain  radius  Rc,  which  is  just  seven  times  the  critical 
radius  of  nucleation  theory.  They  were  able  to  look  at  any 
arbitrarily  shaped  perturbation  by  examining  each  of  its 
harmonics.  If  any  of  these  harmonics  evolve  with  a  faster 
growth  rate  than  the  smooth  interface,  then  this  interface 
will  become  unstable.  The  perturbation  method  of  Mullins  and 
Sekerka  has  proven  to  be  a  most  valuable  tool  in  the  develop¬ 
ment  of  the  theory  of  interface  stability. 

Coriell  and  Parker  [55]  extended  the  results  of  Mullins 
and  Sekerka  for  the  sphere  to  the  case  of  cylinder.  Both 
perturbations  in  circular  cross-section  shape  were  considered 
as  well  as  perturbations  along  the  length  of  the  cylinder.  For 
a  perturbation  of  the  form  r  =  R  +  6  cos(Kcj))  ,  where  k  is  a 
positive  integer,  the  cylinder  was  found  to  be  stable  when 
its  radius  was  less  than  a  critical  radius  R  ,  and  unstable 
when  greater  than  R  ,  analogous  to  the  case  of  a  sphere. 

Tarshis  [56]  followed  the  perturbation  procedure 
originally  proposed  by  Mullins  and  Sekerka  [53]/  but  modified 
it  slightly  to  incorporate  a  coupling  equation  which  provided 
a  means  of  systematic  inclusion  of  other  physical  factors. 

The  results  of  his  study  indicated  that  interface  attachment 
kinetics  can  markedly  affect  the  stability  criterion. 

Tarshish  considered  the  stability  of  a  single  protuberance  and 
not  of  the  individual  harmonics  making  up  some  arbitrary 
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distortion, as  in  the  perturbation  method. 

Kotler  [57]  analyzed  the  stability  of  the  dendrite 
stem,  employing  a  right  circular  cylinder  as  its  physical 
representation.  He  also  treated  the  stability  of  the  para¬ 
boloid  of  revolution  with  the  major  emphasis  on  the  region 
near  the  tip.  In  each  of  these  cases  he  took  into  account 
the  interface  curvature  effects  as  well  as  attachment 
kinetics.  He  concluded  that  the  region  of  greatest  instabil¬ 
ity,  for  most  metals  and  ice,  appeared  to  be  less  than  one 
tip  radius  behind  the  dendrite  tip. 

To  summarize,  it  would  seem  that  at  present  the  theore¬ 
tical  methods  have  arrived  at  a  point  at  which  improvements 
are  very  difficult  to  achieve.  It  may  however  be  that  pro¬ 
gress  can  be  made  using  some  form  of  approximate  numerical 
method.  Since  the  interface  shape  is  always  curved,  and 
boundary  conditions  are  complex,  the  finite  element  method 
seems  to  be  an  appropriate  choice.  The  finite  element  method 
has  been  successfully  used  to  solve  various  partial  differen¬ 
tial  equations  of  fluid-mechanics  and  heat  transfer,  and 
fundamentally  the  dendrite  problem  seems  to  be  similar  in 
form.  The  next  chapter  describes  the  proposed  method  in 
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CHAPTER  3 


FUNDAMENTAL  THEORY 

In  this  chapter  the  fundamental  theory  on  which 
this  research  is  based  will  be  presented.  The  first  portion 
deals  with  the  aspects  of  thermodynamics  that  have  special 
relevance  to  the  crystallization  process,  while  the  later 
portions  deal  with  the  mathematical  formulation  of  the 
problem.  A  brief  description  of  the  finite  element  method 
is  also  presented  for  the  benefit  of  those  unfamiliar  with 
the  method. 

3 . 1  The  Physics  of  Freezing 

In  most  engineering  literature  the  process  by  which 
a  liquid  solidifies  is  described  as  occurring  in  three 
heat  loss  stages.  First  the  liquid  which  was  initially  at 
some  temperature  above  freezing  cools  to  its  equilibrium 
freezing  temperature.  Then  with  further  heat  loss,  solidifi¬ 
cation  begins  and  continues  for  some  time  (depending  on  the 
latent  heat  of  fusion  and  the  rate  of  heat  loss)  at  constant 
temperature.  Finally,  when  all  the  liquid  has  turned  to 
solid,  the  solid  begins  to  cool.  For  example,  suppose  1  gm 
of  water,  initially  at  40°C  and  atmospheric  pressure  is  to 
be  cooled  at  a  rate  of  heat  loss,  say  1  cal/min.  The  cooling 
process  can  be  depicted  ideally  by  Figure  3.1.  The  cooling 
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continues  at  a  constant  rate  for  40  minutes,  at  the  end  of 
which  the  temperature  is  0°C.  With  further  heat  loss,  ice 
begins  to  form  and  continues  forming  for  80  minutes  (the 
heat  of  fusion  of  ice  is,  80  cal/gm) ,  during  which  time  the 
temperature  remains  at  0°C.  At  the  end  of  the  80  minutes, 
cooling  proceeds  once  more,  but  now  at  a  rate  of  approximate¬ 
ly  2  °C/min  because  the  specific  heat  of  ice  is  approximately 
0.5  cal/gm°C. 

In  reality,  when  a  quantity  of  liquid  is  cooled,  freez¬ 
ing  does  not  occur  in  the  ideal  way  shown  in  Figure  3.1. 
Actually,  if  there  is  no  solid  phase  present  a  liquid  can  be 
cooled  substantially  below  its  equilibrium  fusion  temperature 
before  freezing  occurs.  Such  a  liquid  is  in  a  metastable 
supercooled  state.  It  is  found  experimentally  that  the 
maximum  supercooling  attainable  increases  with  increasing 
purity  of  the  liquid,  particularly  as  far  as  foreign  particles 
in  suspensions  are  concerned,  and  that  it  is  easier  to  super¬ 
cool  small  droplets  than  larger  volumes  of  liquid.  Both 
these  effects  are  easily  accounted  for  by  the  assumption 
that  foreign  solid  particles  assist  the  freezing  transition, 
division  of  a  sample  of  liquid  into  many  small  droplets 
effectively  isolating  the  most  active  foreign  particles  into 
a  small  fraction  of  the  total  number  of  droplets.  For  example, 
we  may  easily  cool  tap  water  to  -4  to  -7°C,  while  water  drop¬ 
lets  can  exist  down  to  -40°C  [58].  Except  in  the  case  of 
materials  which  form  glasses,  it  is  not  possible  to  maintain 
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Time  (minutes ) 


Figure  3.1 


Ideal  Cooling  Curve  for  Water 
(From  Knight  [58] ) 
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the  supercooled  state  indefinitely  or  to  achieve  more  than  a 
limited  degree  of  supercooling  before  spontaneous  freezing 
occurs.  Figure  3.2  shows  a  more  realistic  graph  of  cooling 
process  of  1  gm  of  water  at  1  atmospheric  pressure  when 
cooled  at  a  rate  of  1  cal/minute. 

This  behaviour  can  be  very  simply  understood  from 
consideration  of  the  process  of  freezing.  From  Figure  3.3, 
which  shows  the  free  energy-temperature  diagram  with  a  P-T 
diagram,  it  is  clear  that  at  any  temperature  below  that  of 
equilibrium,  the  free  energy  of  the  solid  is  less  than  that  of 
the  liquid,  thus  it  is  energetically  favourable  for  the 
liquid  to  change  to  the  crystalline  state.  This  can  not  be 
accomplished  discontinuously ;  first  a  very  small  volume  of 
liquid  must  crystallize  and  this  crystal  must  then  grow  until 
all  the  liquid  has  frozen.  A  small  crystal  embryo  is,  however, 
in  an  energetically  unfavourable  state  because  of  its  very 
large  surf ace-to-volume  ratio  and  the  positive  free  energy 
associated  with  its  interface  with  the  liquid.  There  is 
thus  a  free  energy  barrier  to  be  overcome  before  freezing  can 
commence  and  this  barrier  can  only  be  surmounted  by  a  nuclea- 
tion  process  depending  upon  thermal  fluctuations.  The  stability 
of  a  small  crystal  embryo  may  be  enhanced  if  it  grows  closely 
upon  an  insoluble  foreign  particle  and  in  this  case  the  nucle- 
ation  process  is  termed  heterogeneous.  If  there  are  no  foreign 
particles  or  surfaces  present,  then  freezing  must  commence 
by  a  process  of  homogeneous  nucleation  within  the  pure  liquid 
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Figure  3.2 


A  Less  Ideal  Cooling  Curve 
for  Water  Showing  Supercooling 
(From  Knight  [58] ) 
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igure  3.3  Free  Energy  -  Temperature  Diagram 

with  a  P-T  Diagram ,  Showing  How 
the  Latter  is  Related  to  Former 
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itself.  This  study  deals  only  with  the  pure  liquids,  thus 
only  the  homogeneous  nucleation  is  of  concern  here. 

From  thermodynamic  treatment  of  equilibrium  across 
a  curved  interface,  it  can  be  shown  that  in  a  pure  supercooled 
melt  a  small  crystal  will  melt  while  a  large  crystal  can 
continue  to  exist  and  grow.  The  proper  criterion  for 
"smallness",  in  this  sense,  is  the  radius  of  curvature.  If 
this  is  less  than  some  value,  the  critical  radius,  that 
depends  on  the  temperature,  melting  occurs;  if  larger,  it 
does  not.  The  critical  radius  can  be  expressed  as:  [3] 


* 


r 


2YTf 

LAT 


(3-1) 


where  Tf  is  the  equilibrium  temperature,  Y  is  the  surface 
free  energy  per  unit  area,  L  is  the  latent  heat  of  fusion  and 
AT  is  the  amount  of  supercooling. 

The  basic  process  to  be  considered  after  nucleation 
is  the  process  of  growth.  The  rate  at  which  freezing  pro¬ 
ceeds  is  controlled  by  the  rate  of  removal  of  the  latent 
heat.  If  it  is  not  removed,  the  temperature  rises  to  the 
point  at  which  no  more  freezing  can  occur  that  is ,  the 
freezing  point. 

The  latent  heat  is  removed  by  conduction;  this  may  be 
through  the  crystal  (as  in  Figure  3.4(a))  or  into  the  liquid, 
if  this  is  supercooled  and  therefore  at  a  lower  temperature 
than  the  interface  where  freezing  occurs  (Figure  3.4(b)) . 

This  study  is  concerned  only  with  the  latter  case. 
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(b)  Heat  Removal  Through  the  Liquid 
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When  the  liquid  is  supercooled  and  freezing  proceeds 
by  rejecting  latent  heat  into  the  liquid,  the  shape  of 
interface  is  not  smooth.  This  can  be  explained  as  follows. 
Figure  3.5  shows  the  schematic  of  removal  of  latent  heat 
through  the  liquid  when  freezing  proceeds  in  a  supercooled 
melt.  On  a  flat  interface  r  =  <=,  the  latent  heat  dissipates 
into  a  prism  of  liquid  [Figure  3.5(a)],  while  the  heat  from 
a  curved  interface  dissipates  into  a  certain  solid  angle 
[Figure  3.5(b)]  and  hence  dissipates  at  a  faster  rate.  Thus, 
any  part  of  the  crystal  interface  which  gets  ahead  of  neigh¬ 
bouring  regions  is  in  a  more  favourable  position  to  dissipate 
latent  heat  of  crystallization  and  so  its  growth  rate  is 
enhanced.  This  situation  leads  to  instability  and  dendritic 
crystal  growth  with  dendrite  arms  branching  out  into  the  liquid. 

3 . 2  Physical  Model  and  Assumptions 

To  begin  simply,  this  study  is  concerned  with  the 
initial  growth  development  of  a  solid  phase  in  a  pure  super¬ 
cooled  melt.  The  geometries  of  interest  are  the  following: 

(a)  Circular  cylinder 

(b)  Parabolic  cylinder 

(c)  Paraboloid  of  revolution 

(d)  Spheroidal 

Particular  interest  in  the  shapes  indicated  arises 
because  considerable  experimental  data  and  analytical  results 
are  available  on  growth  rates  and  stability  of  these  shapes 
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for  growth  in  supercooled  liquids. 

It  is  assumed  that  the  rate  controlling  process  is 
thermal  diffusion,  and  that  interface  kinetics  may  be  ignored. 
The  thermal  profile  in  front  of  the  solid  is  illustrated  in 
Figure  3.6.  The  equilibrium  melting  point  of  the  tip  (inter¬ 
face)  is  depressed  an  amount  AT_^  owing  to  the  radius  of 
curvature  of  the  non-flat  interface.  The  amount  of  freezing 
point  depression  due  to  curvature  effect  (Gibbs-Thompson 
effect)  is  given  by 


AT 


r 


(3-2) 


in  which  y  is  the  surface  free  energy,  K  is  the  mean  curvature, 
L  is  the  latent  heat  generated  per  unit  volume  of  solid 
produced  and  T^  is  the  equilibrium  fusion  temperature.  As 
the  solid  grows,  heat  is  released,  which  is  dissipated  into 
the  supercooled  melt.  This  heat  diffuses  down  a  temperature 
gradient  into  the  melt,  resulting  from  the  temperature 
difference  AT  ;  thus 


AT  =  AT  +  AT  (3-3) 

r  « 

The  other  basic  assumptions  to  be  made  are  the 
following : 

1.  Densities  of  solid  phase  (crystal)  and  the  super¬ 
cooled  liquid  phase  are  equal,  thus,  the  moving  boundary  does 


Thermal  Profile  in  Front  of  a  Dendrite 
Tip  Growing  in  a  Melt  Under  Cooled 
an  Amount  AT  =  Tf  -  T^ 


Figure  3.6 
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not  create  any  convection  in  the  liquid  as  the  growth  proceeds. 

2 .  Physical  properties  are  constant  throughout  each 

phase . 


3.  The  temperature  far  from  the  interface  is  uniform 

<TJ  • 

4.  The  solid  front  moves  with  an  average,  constant 
growth  velocity  v.  This  implies  a  quasi-steady  state  process. 
The  quasi-steady  state  implies  that  the  quantities  of  inter¬ 
est,  for  example  temperatures  in  this  case,  do  not  change 
with  time  with  respect  to  a  coordinate  system  moving  with 
velocity  v. 


3.2.1  Discussion  of  Assumptions 

The  assumption  that  conduction  is  the  main  mechanism 
of  heat  dissipation  is  most  important.  In  the  fusion  problem 
convection  can  arise  from  one  of  two  sources.  Firstly,  a 
convection  due  to  displacement  of  the  liquid  occurs  because 
of  the  different  densities  of  the  liquid  and  solid,  and  pg 
respectively.  As  the  freezing  proceeds,  it  can  be  shown  using 
the  conservation  of  mass  that  the  liquid  is  displaced  ahead 
of  the  crystal  with  a  velocity  given  by 


U  =  -e  v  (3-4) 

n  n 

where  v  is  the  normal  velocity  of  the  interface  and 
n 

e  =  (p  -  p  )/p  .  Chambre  [53]  has  shown  that  unless  e  is 

S  X/  X/ 
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appreciable,  with  respect  to  unity,  the  temperature  distri¬ 
bution  in  the  liquid  and  solid  is  essentially  unchanged  from 
the  e  =  0  case.  Since  e  <  0.1  for  any  material  we  may  wish 
to  consider,  the  effect  of  this  'displacement'  convection 
upon  the  temperature  distribution  during  dendritic  growth 
can  be  justifiably  neglected. 

The  other  type  of  convection,  natural  convection, 
occurs  due  to  density  differences  in  the  melt.  It  is  not  so 
clear  that  this  type  of  convection  can  be  neglected  in  all 
cases.  Gilpin  [60]  has  shown  natural  convection  has  a  signi¬ 
ficant  effect  on  the  experimental  growth  of  dendritic  ice 
for  supercoolings  less  than  about  2°C.  These  experiments  were 
done  in  large  vessels.  Whether  or  not  the  same  conclusion 
would  apply  for  dendrites  growing  in  more  confined  spaces 
is  uncertain.  The  assumption  of  constant  material  properties 
is  justifiable  since  the  temperature  range  of  interest  is  not 
very  large. 

The  assumption  of  uniform  temperature  far  away  from 
the  interface  boundary  was  also  made.  Bolling  and  Tiller  [16] 
have  stated  that  in  practice,  if  the  far  field  boundary  is 
more  than  10  a  /v  away  from  the  solid  interface,  where  a^  is 
the  thermal  diffusivity  of  the  liquid  and  v  is  the  dendrite 
velocity  of  growth,  the  temperature  remains  fairly  constant. 

In  this  study  the  far  field  boundary  is  obtained  by  consider¬ 
ing  various  sizes  of  domains  and  choosing  the  one  which  justifies 
the  assumption  of  nearly  uniform  temperature  far  away  from 


the  interface. 


3 . 3  Mathematical  Formulation  and 
Governing  Equations 
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The  core  problem  is  to  solve  the  quasi-steady  state 
heat  equation  for  the  temperature  fields  inside  and  around 
a  crystal,  of  some  predefined  initial  shape,  growing  into 
a  pure  supercooled  melt. 

The  differential  equation  for  heat  transport  to  be 
solved  for  quasi-steady  state  conditions  in  the  moving 
coordinate  system  fixed  to  the  dendrite  tip  is 
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T  +  — 
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3z 
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(3-5) 


where  a  is  the  thermal  diffusivity,  and  v  is  the  growth 
velocity  of  the  dendrite.  The  velocity  v,  is  of  course,  a 
quantity  that  is  ultimately  to  be  determined.  Setting  the 
right  hand  side  of  the  equation  (3-5)  to  zero  gives  the  quasi¬ 
steady  state  approximation  which  will  be  used  throughout  this 
analysis.  That  is,  it  is  assumed  that  the  diffusion  field 
responds  very  quickly  to  changes  in  the  shape  of  the  moving 
interface.  Also,  since  equation  (3-5)  is  to  be  solved  both 
inside  and  around  the  growing  dendrite,  subscripts  s  and  i 
will  be  used  to  denote  solid  and  liquid  respectively. 

Equation  (3-5)  is  to  be  solved  subject  to  two  boundary 
conditions  at  the  solid-liquid  interface.  First,  there  is  a 
thermodynamic  boundary  condition  due  to  the  Gibbs-Thompson 
effect  that  is  written  in  the  form 
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(3-6) 


in  which  is  the  interface  temperature;  T  is  the  fusion 

temperature  of  the  pure  substance,  K  is  the  curvature  of  the 
surface;  y  is  the  surface  free  energy;  and  L  is  the  latent  heat 
per  unit  volume  the  solid  formed. 

The  second  boundary  condition  at  the  interface  is  the 
heat  balance  condition: 


L  vn  =  <-  k«t  V  Til  +  ks  7  V  '  n  (3-7) 

where  n  is  the  outward  directed  unit  normal  at  a  point  on  the 
surface;  and  v  is  the  'normal  velocity  of  the  surface  at  that 

point.  kp  and  k  are  the  thermal  conductivities  of  the  liquid 

X/  s 

and  solid  phases  respectively. 

In  addition  to  the  interface  conditions,  equals 
(a  constant)  at  a  distance  far  away  from  the  tip  of  the  dendrite. 

3.3.1  Non-Dimens ionaliz at ion 

Equations  (3-5)  through  (3-7)  can  be  non-dimensionali zed 
with  help  of  the  following  substitutions: 

Let 

0  = 


T  -  Tf 
Tf  -  T<* 


r_ 

d 


(3-8a) 


n 


( 3-8b ) 
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(3-8c) 

( 3-8d) 


and 


v  =  v'  v 

c 


K  = 


(3-8e) 


(3-8f) 


where  d  ,  v  and  t  are  characteristic  length,  velocity  and 
c  c  c 

time  respectively.  The  non-dimensional  variables,  n,  £*  v* , 
0  and  K'  are  respectively  length  in  r  direction,  length  in  z 
direction,  velocity,  temperature  and  curvature. 

The  differential  equation  (3-5)  now  becomes 


v2e 


V  '  V 
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36 
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(3-9) 


and  the  interface  conditions  become 


T  y  K' 

9if  =  '  dc  L  ' 


and 
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(3-10) 


(3-11) 


By  choosing  the  characteristic  length  to  be 
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d  =  - = - 

c  L  (T  -  T  )  ' 

r 

equation  (3-10)  reduces  to  0  =  -K* .  Further,  equation 

(3-11)  can  be  simplified  by  letting 


ka  (Tf  -  T~> 
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L  d  v 

c  c 


=  1. 


(3-13) 


This  determines  the  characteristic  time  and  velocity; 


and 
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(3-14) 


(3-15) 


Also,  define  Stefan  number  in  the  liquid 
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and  the  Stefan  number  in  the  solid  as 


(3-16) 
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Then  to  get  a  complete  temperature  distribution  we  need  to 
solve  the  following  differential  equations  for  the  domains 
given  in  Figure  (3-7) : 
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(i)  Liquid  region 
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with  boundary  conditions 


e,  =  -i 


on  S. 
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0  £  -  0if  =  -K'  (n ,£)  on  S 


(ii)  Solid  region 
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with  boundary  conditions 


0  =  0 .  _  on  S  . 

s  if  4 


30 


3n 


=  0 


on  S5,  S6 


and  the  interface  heat  balance  condition  given  by 


0 .  .  =  0„  =  0 

if  £  s 


(3-18) 

(3-19) 

(3-20) 

(3-21) 

(3-22) 

(3-23) 

(3-24) 

(3-25) 


and 
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3 . 4  Variational  Formulation 

Often  boundary  value  problems  have  different  but 
equivalent  formulations — a  differential  formulation  and  a 
variational  formulation.  In  the  first,  a  partial  differential 
equation  is  written,  and  its  direct  solution  is  attempted 
subject  to  given  boundary  conditions.  In  the  second,  the  aim 
is  to  find  the  unknown  function  or  functions  which 

make  stationary  a  functional  or  system  of  functionals 
subject  to  the  same  boundary  conditions.  The  two  problem 
formulations  are  mathematically  equivalent  because  the 
functions  that  satisfy  the  differential  equations  and  their 
boundary  conditions  also  extremize  or  make  stationary  the 
functionals.  Detailed  procedures  for  formulating  many  problems 
of  engineering  and  applied  physics  by  variational  approach 
can  be  found  in  texts  by  Mikhlin  [61] ,  Schecter  [62] ,  and 
Kantorovich  and  Krylov  [63] . 

To  formulate  a  two-dimensional  boundary  value  problem 
of  the  type  given  in  the  last  section  by  a  variational 
approach,  consider  the  following  integral 


1(0) 
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(3-27) 


where  D  is  the  entire  domain  under  consideration.  The  boundary 
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conditions  are  exactly  the  same  as  those  used  in  the  boundary 
value  problem.  The  object  is  to  find  a  function  0(n,£), 
continuous  in  the  region  D  together  with  its  partial  deriva¬ 
tives  of  the  first  and  second  orders,  satisfying  the  specified 
boundary  conditions  and  giving  the  integral  1(0),  Eq.  (3-27) 
a  minimum  value. 

Let  0(n,£)  be  a  function  which  minimizes  the  integral 
1(0) ;  let  us  consider  the  value  of  the  integral  for  the  function 


0(n,o  =  0  (n ,  £)  +  a  0  (n,€) 


(3-28) 


where  0(n,£)  also  satisfies  the  boundary  conditions  and  a  is 


a  parameter.  Since  both  0  and  0  satisfy  the  boundary  conditions, 
3(rif£)  must  vanish  on  the  portion  of  the  boundary  on  which  0 
is  specified. 


To  minimize  1(0),  a  necessary  condition  is  the  vanish¬ 


ing  of  the  first  variation  of  1(0),  which  is  61(0).  The 


requirement  is  therefore 


61(0) 
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(3-29) 


where  1(0  +  a3)  is  given  by 


dnd£ 


(3-30) 
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Differentiating  (3-30)  with  respect  to  a  and  setting  a  =  0 
gives 


61(0) 


ffro  M  1$  ,  o  i§.  !§. 
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(3-31) 


In  the  above  expression  for  61(0),  the  first  two  terms  can 
be  transformed  by  means  of  Green's  formula,  and  equation  (3-31) 
can  be  written  as 


61(0)  =  -2  {  jf  ( V 2 9  +  g)  BdridC  +  J~ [  -  d5  -  |f  dn ]  0  }  =  0 

D  r  dn 

(3-32) 

where  F  is  the  contour  of  the  domain.  As  the  area  integrals 
and  the  line  integrals  are  independent  and  that  $(n,£)  is 
arbitrary,  the  result  is, 


0  +  g  =  0 


on  D 


(3-33) 


and 


(3-34) 


However,  on  the  portion  of  F  (say  F^) ,  on  which  0  is  specified, 
3(n,£)  =  0,  and  on  the  remaining  portion  (say  T  ,  equation 
(3-34)  gives 
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fide 

8n 


96  ^ 

H  dn 


=  0 


on  r 
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(3-35) 


Equation  (3-35)  can  be  reduced  to  the  form 


90 

9n 


tt—  =  0 


on  r 


(3-36) 


when  n  is  the  unit  normal. 

Thus,  we  can  see  that  the  minimization  of  the  integral 
represented  by  equation  (3-27)  results  in  recovering  the  boundary 
value  problems  arrived  at  in  section  3.3  provided  an  approxima¬ 
tion  is  made  such  that 


g 


£ 


90 


v '  Ste 


£ 


£  9£ 


for  the  liquid  region  and 


v '  Ste 

s 


(3-37) 


(3-38) 


for  the  solid  region.  It  is  seen  that  the  functiona.l  does  not 
lead  to  the  partial  differential  equation  if  g  is  taken  as  a 
function  of  0.  However,  the  approximation  considering  g  as 
only  a  function  of  position  is  justified  because  the  terms 
described  by  Eq.  (3-37)  and  (3-38)  are  introduced  iteratively 
in  the  computer  program. 

A  similar  analysis  can  be  made  for  the  case  where 
the  growing  dendrite  has  an  axis  of  symmetry.  The  integral 
to  be  minimized  for  arriving  at  the  boundary  value  problems 
represented  by  equations  (3-18)  and  (3-22)  for  the  axisymmetric 


case  is 
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T/n>  _  fTr,9e,2  .302 

1  ( 0 )  “  JJ  [  (3^)  +  (j^)  -  2  g  0]n  drid^  (3-39) 

where  the  growing  dendrite  is  symmetric  about  the  £  axis. 

3 . 5  The  Finite  Element  Method 

The  finite  element  method  is  one  of  the  newest  and 
most  popular  numerical  techniques  for  solving  the  differential 
equations  of  engineering  and  applied  physics.  This  method  has 
been  used  extensively  in  recent  years  because  it  has,  in 
general,  several  outstanding  advantages.  Some  of  the  main 
ones  include: 

1.  Irregularly  shaped  boundaries  can  be  approximated 
using  elements  with  straight  sides  or  matched  exactly  using 
elements  with  curved  boundaries.  The  method,  therefore,  is 
not  limited  to  "nice"  shapes  with  easily  defined  boundaries. 

2.  The  size  of  the  elements  can  be  varied.  This 
property  allows  the  element  grid  to  be  expanded  or  refined  as 
the  need  arises. 

3.  The  material  properties  in  adjacent  elements  do 
not  have  to  be  the  same.  This  allows  the  method  to  be  applied 
to  non-homogeneous  and  anisotropic  configurations  with 
relative  ease. 

4.  Once  a  computer  program  has  been  developed,  all 
problems  for  which  the  equations  are  the  same  can  be  solved 


•I 
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simply  by  supplying  the  computer  with  appropriate  boundary 
coordinates  and  material  properties. 

The  finite  element  method,  when  applied  to  heat 
transfer  problems,  generally  consists  of  the  following  steps: 

1.  Discretization  of  the  domain  under  study;  defining 
the  nodal  points  and  elements. 

2.  Defining  the  element  function  for  a  single  element 
and  evaluation  of  the  matrices  of  the  elements. 

3 .  Assembling  the  complete  matrix  of  the  continuum 
and  application  of  the  boundary  conditions. 

4.  Solution  of  the  resulting  system  of  equations  to 
obtain  nodal  temperatures. 

5.  Calculation  of  any  other  functions  (temperature 
gradients)  based  on  the  nodal  temperatures. 

The  first  step  of  discretizing  the  domain  involves 
the  decision  as  to  the  number,  size,  and  shape  of  the  elements 
used  to  model  the  domain.  This  decision  has  to  be  based  on 
physical  reasoning  and  experience.  For  example,  the  element 
size  can  be  decreased  in  the  areas  where  the  desired  result 
may  vary  quite  rapidly  (high  gradient  values) ,  and  can  be 
increased  in  regions  where  the  desired  result  is  relatively 
constant. 

Each  element  is  then  analyzed  separately  and  its  prop¬ 
erties  are  generally  derived  from  the  minimization  of  the 
functional  or  Galerkin  type  expression  governing  the  problem, 
after  choosing  a  set  of  functions  to  define  uniquely  the 
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temperature  within  each  "finite  element".  These  approximate 
functions  have  to  satisfy  the  admissibility  and  completeness 
conditions  for  the  problem.  Admissibility  implies  continuity 
of  essential  variables  between  elements,  and  that  the  order 
of  the  expansion  is  such  that  the  terms  are  well  defined  in 
the  variational  statements;  for  completeness,  when  the  elements 
tend  to  be  infinitely  small  and  hence  derivatives  inside  the 
variational  statement  tend  to  be  constant,  the  approximate 
function  must  represent  this  constant  derivative  condition. 

If  these  conditions  are  satisfied  the  solution  will  converge 
to  its  correct  value  as  the  total  number  of  elements  is  in¬ 
creased  . 

3.5.1  The  Finite  Element  Method  Applied  to 
the  Dendritic  Growth  Problem 

In  variational  form  the  temperature  distribution 

problem  to  be  solved  is  that  of  minimizing  the  functional 
given  by  equation  (3-27)  for  the  two-dimensional  problem  and 

for  the  axisymmetric  case  the  functional  given  by  equation 
(3-39)  is  to  be  minimized.  Since  the  procedures  to  be  followed 
in  the  finite  element  formulation  are  the  same  for  two-dimen¬ 
sional  and  axisymmetric  problems,  only  the  axisymmetric  case 
is  described  here  in  detail.  The  main  difference  comes  in 
evaluating  the  element  matrices,  which  are  given  in  Appendix  2. 

Consider  the  case  of  an  axisymmetric  dendrite  growing 
in  a  supercooled  liquid.  The  ultimate  objective  is  to  deter¬ 
mine  the  temperature  distribution  inside  and  around  the 


dendrite,  and  the  temperature  gradients  at  the  interface  to 
calculate  the  growth  velocity  of  the  dendrite  under  given 
conditions . 

First  the  region  under  consideration  is  divided  into  ' 
a  number  of  elements.  The  type  of  the  element  chosen  for  this 
study  is  a  triangular  element  with  six  nodes,  as  shown  in 
Figure  (3.8).  In  order  that  the  element  mesh  may  be  generated 
automatically  by  the  computer  on  supplying  a  minimum  amount 
of  input  data,  the  elements  are  chosen  in  a  block  pattern 
as  shown  in  Figure  (3.9).  The  node  numbering  system  also 
follows  a  distinct  pattern.  The  numbering  system  of  nodes 
and  the  pattern  for  elements  allowed  for  easy  changes  in  total 
number  of  nodes  and  the  size  of  individual  elements.  Details 
of  the  automatic  mesh  generation  procedure  are  given  in 
Appendix  3 . 

The  six  node  triangular  element  enables  a  quadratic 
variation  of  the  temperature  within  the  element.  The  quadra¬ 
tic  element  is  known  to  give  good  results  with  very  little 
extra  computational  work  compared  to  the  linear  element,  and 
without  the  programming  complications  present  in  other  higher 
order  elements  or  the  isoparametric  elements. 

The  generation  of  the  interpolation  function  for  a  tri¬ 
angular  element  is  simplified  considerably  if  one  works  with  the 
area  coordinates.  The  area  coordinates  are  defined  by 


i  =  1/2,3 


(3-40) 
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Figure  3.8  The  Six-node  triangular  Element 

with  area  co-ordinates  of  the  Nodes 
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A  Node  Pattern  for  Automatic  Node 
Numbering 


Figure  3.9 
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where  A  is  the  area  of  the  triangular  element  and  are 
the  areas  of  the  subtriangles  as  shown  in  Figure  (3-8). 
Thus ,  the  side  connecting  nodes  1  and  2  is  described  by 
ip 3  =  0  etc.  We  have  then  three  area  coordinates  (ip  ) 

but  only  two  are  independent  since 


^1+^2  +  ^3  =  1 


(3-41) 


the  mth 

0?  (i  = 


The  quadratic  variation  of  the  temperature,  0m,  in 
element,  expressed  in  terms  of  the  six  nodal  values 
1  to  6)  is 


0 


m 


(3-42) 


where 


X±  =  <  ip1(2ip1  -  1),  ^2(2^2  *'1)'  ^3  ( 2l^3  ‘  1}'  4  *1*2' 


4  ^2^3  '  4^3^1  > 


(3-43) 


Since  r)  is  a  linear  function  in  the  r|-£  plane,  it  is  given 
exactly  by 


n 


(3-44) 
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on  any  triangle  in  the  r|-£  plane,  provided  ri^  are  the  vertex 
radial  coordinates.  Hence,  the  functional  (3-39)  may  be 
expressed  as  the  weighted  sum  of  three  portions 


m  O  ti-i  ^  a  ^ 

1(0)  =  I3  ru  Jf  <1^  [<!£-)  +  (ff-)  -  2g  0m]  dndC 

i=l  1  D  1  3n  35 


(3-45) 

Next,  the  derivatives  with  respect  to  n  and  £  may  be  written 


as 


.m 


3  0  y  6  Q  ^  Ti  /  ■  -I  ,  \ 

—  —  =.L  0.  P.  ( i=l  to  6) 

i=l  ii 


(3-46) 


and 


m 


=.E®  0m  P. 
9£  1=1  i  i 


(i=l  to  6) 


(3-47) 


where 


P. 

l 


- -  {  (4ip  . -1)  b_  ;  ( 4^0-l )  b  ;  (4^  -1)  b  ; 

„  m  l  1  ^  ^  j  J 

2A 


4  (ip2b1  +  4  (i^3b2  +  ^3);  4  (ip^  +  ^3^)} 


(3-48) 


(3-49) 


and 


(3-50) 
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and  p i  is  obtained  by  replacing  the  b's  with  a's  in  the 

expression  for  P_^m. 

Let  the  term  g  be  approximated  by  an  interpolation 
expression 


9  =  z  Xp 

i=l  1  1 


(3-51) 


then  the  equation  (3-45)  can  be  written  as 


i(0m)  =  I3  n.m  ff 
i=l  D 


^  t(E  0  m  P  )2  +  (E  0mp)2-2  E6  g  X  em] 
k=l  k  k  j=l  3  j  n=l  n  n 


dr|d£ 


(3-52) 


To  minimize  (3-52),  we  set 


9 i ( em) 


9  0 
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=  E 
i=l 


m 

ni 


7t-  U  6  m  (P,  P  +  P  P  )  -  2  g  v  v  ] 

J  Yi  .  J  k  J  k  J  n=l  ^n*nAk 

D  3 


dpd£  =  0  (3-53) 


Equation  (3-53)  can  be  simplified  considerably  if  all  qualities 
independent  of  the  triangle  size  and  shape  are  evaluated  once 
for  all.  To  assist  in  doing  so,  define 


SA  ”  -  \ 
iJ  k=1  k 


A  /S 

ff  (PiPJ  +  PiPJ>  dnd5  <3-54> 


D 


^  m 
QAij 


and 


2  S!*±  X±  Xj 


dn  dC 


(3-55) 
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With  this  notation ,  the  minimization  equation  (3-53)  becomes 


SA  9  +  QA  G  =  0 


(3-56) 


where  SA  and  QA  are  6x6  matrices  and  9  and  G  are  6x1  vectors. 

The  functional  minimization  problem  has  thereby  been 
reduced  to  a  simple  matrix  equation,  in  which  the  coefficient 
matrices  SA  and  QA  are  easily  assembled.  The  element  matrices 
SA  and  QA  for  axisymmetric  elements  and  S  and  Q  for  two-dimen¬ 
sional  elements  are  listed  in  Appendix  2.  Also  given  is  a 
sample  integration. 

For  a  dendrite  growing  at  constant  speed,  the  second 
part  of  the  above  matrix  equation  contains  terms  involving  v' 

and  (see  equations  3-37  and  3-38).  Since  v"  and  can 

dt,  dt, 

only  be  obtained  after  the  temperature  field  has  been  determin¬ 
ed,  an  iterative  process  is  used.  That  is,  the  first  computa¬ 
tion  is  done  with  Ste  =  0  so  the  vr  term  drops  out.  The 

dt. 


resulting  temperature  field  is  used  to  calculate  v'  and 


de 

d£ 


c\  ft 

and  the  term  q  =  v'  Ste  —=■  is  introduced  as  an  effective 

dt, 

distributed  heat  generation  in  the  second  iteration.  Typically 
5  or  6  iterations  were  required  to  converge  to  a  self-reprodu- 
cing  result.  It  should  be  noted  that  for  growth  of  a  sphere 
or  a  cylinder,  this  problem  does  not  arise  because  Ste  is 
assumed  to  be  zero. 
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Once  the  individual  matrices  for  all  the  elements 
are  obtained,  essentially  a  bookkeeping  operation  is  required 
to  assemble  the  element  matrices  and  apply  the  boundary 
conditions.  The  method  for  assembling  the  element  matrices 
and  application  of  the  boundary  conditions  can  be  found  in 
standard  texts  [64-66] . 

The  global  matrix  arrived  at  by  the  above  mentioned 
assembly  process  is  always  symmetric  and  banded.  The  band 
width  is  proportional  to  the  largest  difference  between  nodes 
in  the  same  element.  Thus,  considerable  saving  in  storage 
requirements  is  realized  if  only  the  diagonal  and  upper  (or 
lower)  diagnonal  elements  are  stored.  As  shown  in  Figure  (3-9) 
the  numbering  of  the  nodes  is  such  that  the  minimum  bandwidth 


is  realized. 


CHAPTER  4 


STABILITY  OF  PARTICLES  GROWING  IN 
A  SUPERCOOLED  MELT 


4 . 1  Introduction 

The  purpose  of  this  chapter  is  to  examine  the 
applicability  of  the  finite  element  method  to  study  the 
stability  of  a  phase  boundary.  As  mentioned  previously, 
the  stability  is  tested  by  introducing  a  perturbation  in  the 
original  interface  shape  and  determining  whether  this  per¬ 
turbation  will  grow  or  decay. 

Because  some  of  the  heat  flow  problems  arising  out 
of  stability  studies  can  be  solved  exactly  by  analytical 
methods,  an  opportunity  exists  for  checking  computer 
modelling  methods.  The  computer  models  and  their  results 
can  be  compared  at  each  stage  with  analytical  results. 

4 . 2  The  Unperturbed  Growth  of  a  Sphere 

From  a  Slightly  Supercooled  Melt 

Consider  the  case  of  a  solid  sphere  (Figure  4.1) 
growing  in  an  originally  uniformly  supercooled  melt.  If  the 
growth  is  assumed  to  be  diffusion  (heat  transfer)  controlled, 
we  need  to  solve  the  time  dependent  diffusion  equation. 
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Figure  4.1 


The  Geometry  of  a  Sphere  Growing  in  a 
Supercooled  Liquid 
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(4.1) 


subject  to  the  boundary  conditions 


Ta  (r  -  «,  t)  =  (4.2a) 

Tg  (r  +  o,  t)  =  finite  (4.2b) 


T£  ( r=R,  t)  =  Ts  (r=R) ,  t)  =  T±f  =  Tf  -  T  fk 


(4.2c) 


and  the  heat  balance  condition 


9T 


_  k  r  — £.] 
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r=R 


9T 

+  k  [ 
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=  L 


r=R 


dR 
dt  ' 


(4.3) 


to  arrive  at  a  formal  solution  for  the  growth  of  a  sphere. 

An  exact  solution  to  the  above  mentioned  time  dependent 

diffusion  (heat  flow)  equation  is  possible  only  if  the  effect 

of  capillarity  (interface  curvature  effect)  is  neglected.  But 

it  is  not  justifiable  in  real  situations  to  neglect  capillarity. 

It  is  reasonable  to  assume  that  if  the  supercooling  is 

small  then  the  temperature  at  any  point  near  the  growing 

9  T 

particle  changes  slowly  with  time  and  the  term 
(4.1)  is  negligible.  Mathematically  speaking ,  if 


in  equation 
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(T  -  T  ) 
_ f _ «_ 

L 


<<  1, 


(4.4) 


then  we  are  justified  in  approximating  equation  (4.1)  by  the 
Laplace  equation  as  was  shown  for  the  case  of  a  growing  dend¬ 
rite  in  Section  3.3. 


' v. 


i 


-  0. 


(4.5) 


The  solution  of  the  Laplace  equation  (4.5)  subject 
to  the  boundary  conditions  (4.2)  gives 

(T  -  T  ) 

T5,  “  t  +  1  r  R-  (4-6) 


we  find  from  equation  (4.6) 
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Z. 


1-^1 

9r 


r=R 


T .  -  T 

if 

R 


(4.7) 


and  the  growth  rate  of  the  sphere  can  be  obtained  from  equation 
(4.3)  to  be 


R 


dR 

dt 


T]  , 


(4.8) 


assuming  T  =  T.  for  all  r  <  R. 

S  1.  L 

The  results  based  on  the  equations  (4. 6-4. 8)  provide 
an  excellent  starting  ground  for  testing  the  accuracy  of  the 
finite  element  method  for  solving  crystal  growth  problems. 
The  finite  element  programme  was  run  for  a  sphere  having  an 
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initial  radius  ranging  from  50  R*  to  250  R* ,  where 
R*  —  2T  T^/  (Tf-T^)  is  the  classical  nucleation  radius  of 
a  sphere.  Since  it  is  not  possible  to  have  the  outer  boundary 
of  the  supercooled  melt  at  infinity  while  using  the  finite 
element  method,  an  approximation  has  to  be  made  that  at  some 
sufficiently  large  finite  distance  away  from  the  interface 
of  the  growing  sphere  the  melt  is  at  a  uniform  temperature,  T  . 
In  the  present  analysis  the  radius  of  the  outer  boundary  was 
taken  to  be  60  times  the  radius  of  the  ice  sphere. 

The  domain  to  be  considered  for  the  finite  element 
solution  of  the  sphere  growth  problem  is  shown  in  Figure 
(4.2).  For  this  it  is  enough  to  consider  just  the  one  quadrant 
as  shown  in  Figure  (4.2). 

Figures  (4. 3-4.5)  and  Tables  (4. 1-4. 3)  present  the 
comparison  of  the  results  obtained  by  the  finite  element 
method  with  the  analytical  solutions.  It  is  evident  that  the 
finite  element  method  approximates  the  distribution  of  temper¬ 
ature  and  temperature  gradient  around  a  growing  sphere  very 
accurately,  and  the  same  is  true  for  predicting  the  growth 
rate  for  spheres  of  different  initial  radii.  Further,  it  is 
seen  that  the, finite  element  method  predicts  the  temperature 
distribution  extremely  accurately  near  the  interface,  while 
the  accuracy  is  not  so  great  far  away  from  the  interface. 

But,  in  stability  problems  and  the  problems  of  dendritic  growth 
we  only  need  the  temperature  distribution  near  the  interface, 
thus  further  justifying  the  assumption  of  considering  an 
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Figure  4.2  The  Domain  to  be  Considered  for  the 

Study  of  the  Growth  of  a  Sphere 
in  Supercooled  Liquid 
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Figure  4.3  A  Comparison  of  the  Exact  Solution 

and  the  Finite  Element  Solution  for  the 
Temperature  Distribution  Around  a 
Sphere  Growing  in  a  Supercooled  Liquid 

(R  =  100  R*;  0 


Along  Horizontal  Axis 
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A  Comparison  of  the  Exact  Solution 
and  the  Finite  Element  Solution  for 
the  Temperature  Gradient  Around  a 
Sphere  Growing  in  Supercooled  Liquid 

(R  =  100  R*) 


Figure  4.4 


Along  the  Vertical  Axis 
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A  Comparison  of  the  Exact  Solution 
and  the  Finite  Element  Solution  for 
the  Temperature  Gradient  Around  a 
Sphere  Growing  in  a  Supercooled  Liquid 

(R  =  100  R*) 


Figure  4.5 


TABLE  4 . 1 


A  COMPARISON  OF  THE  EXACT  SOLUTION  AND  THE  FINITE  ELEMENT 

SOLUTION  FOR  THE  GROWTH  RATE  OF  VARIOUS  POINTS 

ALONG  THE  INTERFACE  OF  A  SPHERE  GROWING  IN 

SUPERCOOLED  LIQUID  (R  =  50  R* ,  k  =0) 

s 


Position  Along  the 
Interface  (Angle  in 
Degrees  from  the 

Z  Axis) 

Finite  Element 

Solution  for  the 

Growth  Rate  (Non- 
Dimensionali zed ) 

Exact  Solution 
for  the  Growth 
Rate  (Non-Dimen- 
sionalized) 

0 

Q. 002460 

0.002506 

9 

0.002506 

0.002506 

18 

0.002502 

0.002506 

27 

0.002503 

0.002506 

36 

0.002504 

0.002506 

45 

0.002506 

0.002506 

54 

0.002507 

0.002506 

63 

0.002508 

0.002506 

72 

0.002506 

0.002506 

81 

0.002507 

0.002506 

90 

0.002506 

0.002506 
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TABLE  4 . 2 

A  COMPARISON  OF  THE  EXACT  SOLUTION  AND  THE  FINITE  ELEMENT 

SOLUTION  FOR  THE  GROWTH  RATE  OF  VARIOUS  POINTS 

ALONG  THE  INTERFACE  OF  A  SPHERE  GROWING  IN  A 

SUPERCOOLED  LIQUID  (R  =  100  R* ,  k  =0) 

s 


Position  Along  the 
Interface  (Angle  in 
Degrees  frcm  the 
Z  Axis) 


Finite  Element 
Solution  for  the 
Growth  Rate  (Non- 
Dimens  ionalized ) 


Exact  Solution 
for  the  Growth 
Rate  (Non-Dimen- 
sionalized) 


0 

0.001245 

0.001253 

9 

0.001261 

0.001253 

18 

0.001262 

0.001253 

27 

0.001264  ; 

0.001253 

36 

0.001264 

0.001253 

45 

0.001264 

0.001253 

54 

0.001265 

0.001253 

63 

0.001264 

0.001253 

72 

0.001263 

0.001253 

81 

0.001263 

0.001253 

90 

0.001263 

0.001253 
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TABLE  4.3 

A  COMPARISON  OF  THE  EXACT  SOLUTION  AND  THE  FINITE  ELEMENT 

SOLUTION  FOR  THE  GROWTH  RATE  OF  VARIOUS  POINTS 

ALONG  THE  INTERFACE  OF  A  SPHERE  GROWING  IN  A 

k 

SUPERCOOLED  LIQUID  (R  =  250*,  ^  =  4 ) 

l 


Position  Along  the 
Interface  (Angle  in 
Degrees  from  the 
Z  Axis) 


Finite  Element 
Solution  for  the 
Growth  Rate  (Non- 
Dimensionalized ) 


Exact  Solution 
for  the  Growth 
Rate  (Non-Dimen- 
sionalized) 


0 

0.000495 

0.000499 

9 

0.000502 

0.000499 

18 

0.000502 

0.000499 

27 

0.000503 

0.000499 

36 

0.000503 

0.000499 

45 

0.000503 

0.000499 

63 

0.000503 

0.000499 

72 

0.000502 

0.000499 

81 

0.000502 

0.000499 

90 

0.000502 

0.000499 

77 


arbitrary  distance  far  away  from  the  interface  to  be  infinitely 
away  from  the  origin.  A  node  that  caused  consistent  problems 
is  the  one  on  the  ice  surface  at  zero  degrees.  The  growth 
rate  calculated  for  this  node  is  consistently  lower  than  that 
for  other  nodes  around  the  sphere.  This  anamoly  is  an  artifact 
produced  by  the  fact  that  this  node  is  at  a  corner  of  the 
solution  domain  where  a  constant  temperature  surface  meets  an 
insulated  surface  (i.e.  a  singularity  exists).  Methods  to  re¬ 
duce  the  error  resulting  from  such  singularities  are  available 
in  the  finite  element  literature  (67,  68,  69]. 


4 . 3  The  Stability  of  a  Sphere  Growing  in 
a  Slightly  Supercooled  Liquid 

4.3.1  Theoretical  Solutions 

We  consider  a  sphere  of  radius  R  slightly  perturbed  by 
a  spherical  harmonic  The  position  of  the  perturbed 

surface  is  given  by 


r  =  R  +  6  Y^m  (0,<f>)  ,  (4.9) 

,6 

where  6  is  a  small  first  order  quantity  (^-  <<  1)  ;  and 
Y n  (0,cj>)  is  a  spherical  harmonic  which  determines  the  shape 
contours  position  along  the  surface  of  the  sphere  as  shown 
in  Table  4.4  and  Figure  4.6.  For  a  small  perturbation,  the 
curvature  of  a  slightly  perturbed  sphere  may  be  approximated  by 


2  Y 


&m 


6 


(Z2  +  Z) 


6 


K 


2 

R 


+ 


(4.1Q) 
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TABLE  4.4 

SOME  SPHERICAL  HARMONICS ,  Yn 

)6ra 

Y1q  =  Cos  0 

Y  ^  =  Sin  0  Sin  0  or  Sin  0  Cos  $ 

Y20  =  0.5  (3  Cos* 2  0  -1) 

Y21  =  3  CoS  0  Sin  0  C°S  ^  or  3  Cos  0  Sin  0  Sin  ^ 

2  2 

Y  9  =  3  Sin  0  Cos  (2  0)  or  3  Sin  0  Sin  (2  0) 

Y30  =  0.5  (5  Cos3 *  0  - 


3  Cos  0) 


A  Sphere  Perturbed  by  the 
Harmonic  Y^g.  The  Figure  is 

Cylindrically  Symmetric  About 
the  Vertical  Axis 


Figure  4 . 6 
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Further,  since  6  is  small  compared  to  R  and  in  linear  stability 
theory  each  harmonic  develops  independently  of  the  others,  the 
velocity  of  the  growth  of  the  perturbed  sphere  is  also  given 
by 


v 


R  +6 


(4.11) 


It  can  be  shown  that  the  temperature  distribution 
in  the  liquid  and  solid  around  a  perturbed  sphere  discussed 
above  is  given  by  [53] 


T£  (r,0,<|>) 


(Tf-TJR  -  2Tfr  {  (Tf-Ta)  R 


Ti  R^  1  £(£+1)}  6Y 
f  £m 

£+1 

r 


+  T 

oc 


(4.12) 


and 


T 

s 


(r ,  dr(p ) 


2I\ 

“  V1-  R-> 


T  (Z+2 )  (il-l)r  6Yfm 

£  +  2 
r)6  +  z 


(4.13) 


Knowing  the  temperature  distribution  from  equations  (4.12)  and 
(4.13),  the  velocity  v  of  each  element  of  the  interface  may  be 
calculated  from 


v 


k 

s 

L 


9  T 

( - — ) 

1  9r 


i .  f . 


8T£ 


f 

i .  f . 


(4.14) 


where  the  subscript  i.f.  indicates  that  the  quantity  is  to  be 
evaluated  at  the  interface.  It  is  assumed  in  the  equation 
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(4.14)  that  the  normal  gradients  on  the  interface  are  the 
same  as  the  radial  gradients. 

The  quantity  6/6  can  now  be  evaluated  using  the 
equation  (4.11)  and  it  can  be  shown  to  be 

6  (A"1)  Tf  K£  Tf~Tcc  r  k 

J  =  - 2 {  ~T - R  [  U+2)  U+l)  +  2+1  (£  +  2)  } 

R  L  f  k£ 

(4.15) 

The  two  criteria  for  stability  are  defined  as 

•  • 

(6/6)/(R/R)  =  0,  (4.16) 

for  absolute  stability  (no  growth  of  perturbation) ,  and 

•  • 

(6/6)/(R/R)  <  1,  (4.17) 

for  relative  stability  (perturbation  grows  no  faster  than 

the  sphere  itself) .  It  is  clear  that  the  quantity 

•  • 

(6/5)/(R/R)  can  be  easily  calculated  from  the  equations 

(4.15)  and  (4.8)  . 

Several  computer  runs  were  made  to  test  the  applica¬ 
bility  of  the  finite  element  method  to  the  problem  of  the 
stability  of  a  slightly  perturbed  sphere  growing  in  an 
originally  uniformly  supercooled  melt.  The  angular  harmonics 
considered  were  Y^q  and  Figures  4.7  through  4.9  and 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  From  the  Z  Axis) 


A  Comparison  of  the  Analytical  Solution 
and  the  Finite  Element  Solution  for  the 
Growth  Rate  of  a  Sphere  Perturbed  by  the 


Harmonic  Y 


30 


(R  =  100  R*; 


IT-  =  D 


Figure  4.7 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  from  the  Z  Axis) 

Figure  4.8  A  Comparison  of  the  Analytical  Solution 

With  the  Finite  Element  Solution  for 
the  Growth  Rate  of  a  Sphere  Perturbed 
by  the  Harmonic  Y^0  (R  =  100  R*; 


0) 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  From  the  Z  Axis) 


Figure  4.9  A  Comparison  of  the  Analytical  Solution 

and  the  Finite  Element  Solution  for  the 
Growth  Rate  of  a  Sphere  Perturbed  by  the 
Harmonic  Y^q  (R  =  100  R*;  ks 

*7 


4) 
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and  Table  4.5  and  4.6  summarize  the  comparison  of  the  finite 
element  calculations  and  the  analytical  results.  It  is  evident 

that  generally  good  agreement  exists  between  the  two  methods, 

•  • 

the  maximum  error  in  the  quantity  6/6  /  R/R  being  about  40 
percent  for  the  case  ks/k&  =  4,  but  less  than  10  percent  when 

ks/k^  =  0.  However,  it  should  be  noted  that  these  errors 

•  • 

numerically  magnified  because  the  calculations  for  6/6  /  R/R 
involve  subtraction  of  quantities  which  are  very  close  to  each 
other.  This  means  that  the  finite  element  results  are  much 
closer  to  the  analytical  results  when  the  heat  transfer  in  the 
solid  side  is  neglected.  This  should  be  expected  because  it  is 
difficult  to  model  the  solid  side  of  the  sphere  very  accurately 
since  the  sizes  involved  are  so  small.  The  accuracy  of  the 
finite  element  calculations  can  be  further  improved  by  increas¬ 
ing  the  number  of  nodes. 

4 . 4  The  Stability  of  a  Solid  Cylinder 
Growing  in  a  Supercooled  Liquid 

The  next  geometrical  shape  considered  is  a  cylinder 
in  an  infinite  liquid.  The  length  of  the  cylinder  is  assumed 
to  be  large  enough  to  avoid  end  corrections. 

4.4.1  The  Unperturbed  Growth  of  a  Cylinder 

Using  arguments  similar  to  those  used  for  the  sphere 
would  result  in  the  Laplace  equation 
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TABLE  4.5 

A  COMPARISON  OF  THE  ANALYTICAL  METHOD  AND  THE  FINITE  ELEMENT 
METHOD  TO  STUDY  THE  STABILITY  OF  A  SLIGHTLY  PERTURBED 
SPHERE  GROWING  IN  A  SUPERCOOLED  LIQUID 
(THE  SPHERE  PERTURBED  BY  THE  HARMONIC  Y^  Q ) 


R/R* 

k  /k 
s'  a 

Analytical 

•  • 

6/6  /  R/R 

Finite  Element 

•  • 

6/6  /  R/R 

25 

0 

0.750 

0.792 

25 

l 

0.583 

0.634 

25 

4 

0.083 

0.116 

50 

0 

0.878 

0.913 

50 

1 

0.796 

0.873 

50 

4 

0.551 

0.695 

100 

0 

0.939 

•0.961 

100 

1 

0.899 

0.942 

100 

4 

0.777 

0.901 

87 


TABLE  4 . 6 

A  COMPARISON  OF  THE  ANALYTICAL  METHOD  AND  THE  FINITE  ELEMENT 
METHOD  TO  STUDY  THE  STABILITY  OF  A  SLIGHTLY  PERTURBED 
SPHERE  GROWING  IN  A  SUPERCOOLED  LIQUID 
(THE  SPHERE  PERTURBED  BY  THE  HARMONIC  Y^) 


R/R* 

k  /k 

s/  % 

Analytical 

Finite  Element 

6/6  /  R/R 

6/6  /  R/R 

25 

0 

i 

1.167 

1.223 

25 

1 

0.542 

0.595 

25 

4 

-1.333 

-1.220 

50 

0 

1.592 

1.766 

50 

1 

1.286 

1.53 

50 

4 

0.367 

0.624 

100 

0 

1.798 

1.914 

100 

1 

1.647 

1.843 

100 

4 

1.191 

1.434 
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7  T*  =  °' 


(4.18) 


subject  to  the  boundary  conditions 


Tf  (r  =  a)  -  T 

36  a 


(4.19) 


TT 

Vr=R)  =  Tf-^i=T.f( 


(4.20) 


and  the  heat  balance  condition 


3T 


i 


3T 


kJ,  [3r  ]  +  ks  [3r  1 

r=R  r=R 


=  L 


dR 

dt 


(4.21) 


The  solution  of  the  above  equation  (4.18)  is 


(r)  =  A  £n(r)  +  B 


(4.22) 


where  A  and  B  are  to  be  evaluated  using  boundary  conditions. 
However,  the  boundary  condition  T  (r=a)  =  T  ,  makes  T  at 

36  QL  36 

r  =  a  infinite,  so  there  is  a  difficulty  in  using  the  Laplace 
equation  with  a  boundary  at  r  =  a .  The  steady  state  cylindrical 
case  can,  therefore,  only  be  solved  in  a  finite  domain.  If 
the  outer  boundary  is  taken  as  R^ ,  the  temperature  distribution 
is 


T,(r) 


T  + 

a 


[1 


£n  (r/R)  , 
£n  (R^/R) J 


(4.23) 


when  equation  (4.23)  is  used  in  the  boundary  condition  for 
the  conservation  of  heat  we  find 
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dR 


rT 


f 


dt 


LR  &n  (R  /R) 


R 


(4.24) 


Figures  4.10  through  4.12  and  Tables  4.7  and  4.8 


present  the  comparison  of  the  finite  element  calculations  and 
analytical  results  for  the  growth  of  a  cylinder  in  a  super¬ 
cooled  liquid.  The  results  show  a  very  close  agreement.  It 
should  be  noted  that  in  the  case  of  a  sphere  the  results 
obtained  from  a  finite  but  sufficiently  large  domain  were 
compared  with  analytical  solutions  based  on  an  infinite  domain 
while  in  the  case  of  the  cylinder  the  results  from  a  finite 
domain  are  compared  with  analytical  solutions  based  on  a  finite 
domain.  Thus,  it  should  be  expected  that  the  finite  element 
results  for  the  cylinder  agree  much  more  closely  with  the 
analytical  results. 

4.4.2  The  Stability  of  a  Perturbed  Cylinder 


Consider  a  cylinder  of  radius  R  slightly  perturbed  by 


a  set  of  spatial  harmonics.  The  position  of  the  perturbed 
surface  is  given  by 


(4.25) 


R  +  6  cos  ( Kef) ) 


r 


where  6  is  a  positive  length  very  small  compared  to  R  and 
K  is  a  positive  integer,  so  that  r  returns  to  its  original 
value  when  4>  changes  by  2 tt  .  Hardy  and  Coriell  [35]  have  shown 
that  to  the  first  order,  if  normal  gradients  on  the  interface 
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A  Comparison  of  the  Analytical  Solution 
and  the  Finite  Element  Solution  for  the 
Temperature  Distribution  Around  a 
Cylinder  Growing  in  a  Supercooled  Liquid 

(R  =  50  R*;  Rx/R  =  60) 


Figure  4.10 


Along  the  Vertical  Axis 
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A  Comparison  of  the  Analytical  Solution 
and  the  Finite  Element  Solution  for  the 
Temperature  Gradient  Around  a  Cylinder 
Growing  in  a  Supercooled  Liquid 

(R  =  50  R*;  Ra/R  =  6°) 


Figure  4.11 
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A  Comparison  of  the  Analytical  Solution 
and  the  Finite  Element  Solution  for  the 
Temperature  Gradient  Around  a  Cylinder 
Growing  in  a  Supercooled  Liquid 

(R  =  50  R*;  R^/R  =  60) 


Figure  4.12 


93 


TABLE  4 . 7 

A  COMPARISON  OF  THE  EXACT  SOLUTION  AND  THE  FINITE  ELEMENT 
SOLUTION  FOR  THE  GROWTH  RATE  OF  VARIOUS  POINTS  ALONG  THE 
INTERFACE  OF  A  CYLINDER  GROWING  IN  A  SUPERCOOLED  LIQUID 
(R  =  1024  R*,  Rx/R  =  6 0 ,  kg/k£  =  4) 


Position  Along  the 
Interface  (Angle  in 
Degrees  frcm  the 
Z  Axis) 


Finite  Elenent 
Solution  for  the 
Growth  Rate  (Non- 
Dintensionalized ) 


Exact  Solution 
for  the  Growth 
Rate  (Non-Dimen- 
sionalized ) 


0 

0.6000 

E  -04 

0.6010 

E  -04 

12 

0.6017 

E  -04 

0.6010 

E  -04 

24 

0.6017 

E  -04 

0.6010 

E  -04 

36 

0.6018 

E  -04 

0.6010 

E  -04 

48 

0.6018 

E  -04 

0.6010 

E  -04 

60 

0.6019 

E  -04 

0.6010 

E  -04 

72 

0.6019 

E  -04 

0.6010 

E  -04 

90 

0.6062 

E  -04 

0.6010 

E  -04 
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TABLE  4.8 

A  COMPARISON  OF  THE  EXACT  SOLUTION  AND  THE  FINITE  ELEMENT 
SOLUTION  FOR  THE  GROWTH  RATE  OF  VARIOUS  POINTS  ALONG  THE 
INTERFACE  OF  A  CYLINDER  GROWING  IN  A  SUPERCOOLED  LIQUID 
(R  =  512  R*,  RA/R  =  60,  k  /k  =  4) 

S  X/ 


Position  Along  the 
Interface  (Angle  in 
Degrees  frcm  the 
Z  Axis) 


Finite  Elenent 
Solution  for  the 
Growth  Rate  (Non- 
Dimensionalized ) 


Exact  Solution 
for  the  Growth 
Rate  (Non-Dimen- 
sionalized) 


0 

0.1199 

E  -03 

0.1202 

E  -03 

12 

0.1204 

E  -03 

0.1202 

E  -03 

24 

0.1204 

E  -03 

0.1202 

E  -03 

36 

0.1204 

E  -03 

0.1202 

E  -03 

48 

0.1204 

E  -03 

0.1202 

E  -03 

60 

0.1204 

E  -03 

0.1202 

E  -03 

72 

0.1204 

E  -03 

0.1202 

E  -03 

90 

0.1212 

E  -03 

0.1202 

E  -03 
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are  considered  the  same  as  radial  gradients 


6 

6 


R 


R 


* 


k  (ir-i)  (i+  j^-) } 


(4.26) 


* 


where  R  is  the  classical  nucleation  radius.  Equations  (4.26) 

and  (4.24)  can  be  used  in  combination  to  obtain  the  quantity 

•  • 

(6/6)/(R/R)  which  determines  the  stability  criterion. 

Several  computer  runs  were  made  using  different  values  of  K, 
R^/R,  and  kg/k^  for  testing  the  stability  of  a  perturbed 
cylinder.  Figures  4.13  through  4.15  and  Tables  4.9  through 
4.10  show  the  comparison  of  the  finite  element  results  with 
the  analytical  results.  Again,  the  results  are  in  good  agree¬ 
ment  and  the  accuracy  is  even  better  when  the  solid  side  heat 
transfer  is  neglected. 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  from  the  Z  Axis) 


Figure  4.13  A  Comparison  of  the  Analytical  Solution 

with  the  Finite  Element  Solution  for  the 
Growth  Rate  of  a  Perturbed  Cylinder 

(R  =  1024  R*; 


K  =  6) 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  from  the  Z  Axis) 


Figure  4.14  A  Comparison  of  the  Analytical  Solution 

with  the  Finite  Element  Solution  for  the 

Growth  Rate  of  a  Perturbed  Cylinder 

k 

(R  =  1024  R  ;  £ —  =  1 /  K  =  6) 

£ 


0.4  0.6 
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Angular  Position  Along  the  Interface 
(Angle  in  Degrees  from  the  Z  Axis) 


Figure  4.15  A  Comparison  of  the  Analytical 

Solution  with  the  Finite  Element 
Solution  for  the  Growth  of  a  Perturbed 
Cylinder  (R  =  512  R*; 
k 
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TABLE  4 . 9 

A  COMPARISON  OF  THE  ANALYTICAL  METHOD  AND  THE  FINITE  ELEMENT 
METHOD  TO  STUDY  THE  STABILITY  OF  A  SLIGHTLY  PERTURBED 
CYLINDER  GROWING  IN  A  SUPERCOOLED  LIQUID 
(THE  CYLINDER  PERTURBED  BY  THE  HARMONIC  K=6) 


R/R*  ks/ka 


Analytical  Finite  Element 

6/6  /  R/R  6/6  /  R/R 


400 

0 

2.845 

3.012 

400 

1 

0.877 

1.013 

400 

.4 

-5.770 

-7.234 

512 

0 

3.317 

3.498 

512 

1 

1.635 

1.785 

512 

4 

-3.413 

o 

00 

• 

1 

1024 

0 

4.160 

4.390 

1024 

4 

0.798 

1.209 

• 

100 


TABLE  4.10 

A  COMPARISON  OF  THE  ANALYTICAL  METHOD  AND  THE  FINITE  ELEMENT 
METHOD  TO  STUDY  THE  STABILITY  OF  A  SLIGHTLY  PERTURBED 
CYLINDER  GROWING  IN  A  SUPERCOOLED  LIQUID 
(THE  CYLINDER  PERTURBED  BY  THE  HARMONIC  K=2) 


R/R* 


VkSL 


Analytical  Finite  Element 

6/6  /  R/R  6/6  /  R/R 


100 

100 

100 

200 

200 

200 

400 

400 


1 

1 

4 

0 

1 

4 

0 

1 


0.752 

0.504 

-0.240 

0.876 

0.753 

0.383 

0.938 

0.876 

0.692 


0.791 

0.612 

-0.430 

0.923 

0.862 

0.501 

0.952 

0.930 

0.783 
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CHAPTER  5 


DENDRITIC  GROWTH 

5 . 1  Introduction 

The  purpose  of  this  chapter  is  to  demonstrate  that 
the  finite  element  method  can  be  used  effectively  to  calculate 
the  velocity  of  growth  of  freely  growing  dendrites  in  super¬ 
cooled  liquids.  The  first  portion  of  this  chapter  deals  with 
the  selection  of  an  appropriate  domain  for  the  finite  element 
study,  while  the  latter  portion  presents  the  results  obtained 
using  this  domain  for  various  systems  (i.e.  dendritic  growth 
in  ice,  tin,  succiononitrile  etc.).  The  results  thus 
obtained  are  compared  with  the  available  analytical  and 
experimental  results.  The  shapes  studied  are  paraboloids  of 
revolution  and  parabolic  cylinders. 

5 . 2  Selection  of  Domain  to  be  Studied 

As  described  previously  in  Chapter  2 ,  most  of  the 
analytical  studies  of  dendritic  growth  consider  a  dendrite 
growing  in  a  supercooled  liquid  infinite  in  extent.  However, 
in  practice  the  domain  is  always  finite.  Also,  since  both 
the  paraboloid  of  revolution  and  parabolic  cylinder  have  an 
axis  of  symmetry,  the  domain  to  be  studied  can  be  represented 
by  Figure  5.1.  The  problem  now  is  to  select  various  dimensions 
in  the  Figure  5.1  so  that  the  finite  element  results  become 
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element  study  of  dendritic  growth 
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as  realistic  as  possible.  To  achieve  this  goal  the  following 
criteria  should  be  met:  (1)  the  outer  boundary  should  be 

close  to  an  isotherm  with  temperature  equal  to  T  ,  (2)  the 

0 1 

domain  should  be  large  enough  so  that  the  calculated  inter¬ 
face  velocity  changes  very  little  by  increasing 
domain  any  further,  (3)  the  number  of  nodes  should  be 
optimized  in  a  way  that  the  results  are  reasonably  accurate 
yet  the  computer  time  required  is  not  excessive. 

To  avoid  a  costly  and  long  trial  and  error  procedure 
to  obtain  an  appropriate  domain  some  guidelines  were  adopted 
from  the  available  literature.  Bolling  and  Tiller  [16]  have 
shown  that  a  melt  may  be  considered  infinite  with  respect  to 
the  dendrite  when  its  boundaries  are  'v  10  cx^/v  away  from  the 
body  of  the  growing  dendrite  and  it  is  sufficiently  long  for 
the  steady  state  to  be  reached.  Kotler  and  Tarshis  [39]  have 
noted  that  a  paraboloid  of  revolution  is  an  excellent  approxi¬ 
mation  for  the  shape  over  a  distance  on  the  order  of  10R  back 
from  the  tip.  Holzman  [46]  has  presented  a  contour  map  of 
isotherms  in  and  around  a  growing  dendrite  which  can  be  used 
to  define  a  domain  to  be  studied.  An  approximate  value  of 
the  tip  curvature  of  a  dendrite  for  various  systems  can  be 
easily  determined  from  the  available  literature. 

Keeping  the  above  mentioned  guidelines  in  mind,  an 
initial  domain  was  selected.  The  finite  element  method  as 
described  in  Chapter  3  was  used  to  obtain  the  temperature 
field  inside  and  around  a  growing  dendrite.  The  assumptions 
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and  mathematical  formulation  were  the  same  as  those  described 
in  Chapter  3.  Sensitivity  studies  were  made  to  select  the 
most  appropriate  domain  to  be  studied.  The  finite  element 
grid  was  such  that  the  mesh  was  finer  near  the  dendrite  tip 
(where  temperature  gradients  are  large)  and  coarse  away  from 
the  tip  (where  temperature  gradients  are  small).  Figure  5.2 
presents  the  effect  of  the  size  of  the  dimension  on  the 
tip  velocity  of  the  dendrite.  It  is  clear  that  when  the 
domain  is  small  (i.e.  when  the  boundary  with  the  temperature 
T  is  not  very  far  from  the  tip)  the  driving  force  is  such 

LX 

that  the  dendrite  will  grow  with  a  high  velocity.  However, 
as  the  domain  is  made  larger  the  tip  velocity  tends  to  reach 
a  limiting  value.  To  fix  the  dimension  it  can  be  argued  on 
physical  grounds  that  R^  should  be  larger  than  R2  so  that  the 
temperature  gradients  near  the  outer  boundary  are  nearly  the 
same.  After  only  a  few  trials  it  was  found  that  to  study  the 
growth  of  dendrites  with  the  shape  of  paraboloids  of  revolution 
the  most  appropriate  domain  was  such  that  R^  =  10R,  R2  =  20R, 
and  R^  =  4 OR.  The  temperature  field  obtained  from  this  domain 
further  justified  its  appropriateness.  Figures  5.3  and  5.4 
present  the  temperature  field  along  the  vertical  and  horizontal 
axes  respectively.  It  is  clear  that  the  temperature  gradients 
near  the  outer  boundary  in  both  figures  are  very  small  and  equal 
(ideally  the  temperature  gradient  at  the  outer  boundary  should 
approach  to  zero);  while  the  gradients  near  the  interface  are 
such  that  the  tip  velocity  is  greater  than  the  velocity  at  the 
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Figure  5. 2 


Plot  to  determine  the  effect  of  the  finite 
element  danain  size  on  the  growth  rate  of 
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Non-Dimensional  Distance  from  the  Origin 
Along  the  r  Axis 


Figure  5.3  Plot  of  the  temperature  distribution 

around  a  dendrite  growing  in  a  super¬ 
cooled  liquid 
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Distance  from  the  Tip  Along  Z  Axis 


Figure  5.4  Plot  of  the  Temperature  Distribution 

around  a  Dendrite  Growing  in  a 
Supercooled  Liquid 
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bottom  end  of  the  dendrite.  This  is  very  true  in  real 
situations.  Each  side  of  the  domain  (liquid  and  solid)  was 
divided  into  200  triangular  elements.  It  was  found  that  by 
increasing  the  number  of  nodes  any  further  the  accuracy  did 
not  change  significantly. 

5 . 3  Dendritic  Growth  Velocity 
5.3.1  Isothermal  Dendrite 

Ivantsov's  calculation  has  already  been  mentioned  in 
Chapter  2.  It  was  easily  possible  to  simulate  his  model  by 
the  finite  element  method  since  the  only  change  needed  in  the 
problem  formulation  was  to  eliminate  the  Gibbs-Thompson  effect. 
The  main  characteristic  of  the  Ivantsov  model  is  that  for 
different  tip  curvatures  chosen  the  tip  velocity  is  such  that 
the  product  vR  is  constant.  Also,  the  axial  component  of 
surface  growth  remains  constant  since  the  dendrite  retains 
its  shape  as  it  grows  according  to  the  Ivantsov  model.  Figure 
5.5  presents  the  comparison  of  the  tip  velocity  obtained  by 
the  finite  element  method  and  the  analytical  method  for  an 
isothermal  dendrite.  The  results  are  in  very  close  agreement. 
Also,  it  is  seen  that  the  product  vR  remains  almost  constant 
in  the  finite  element  calculations.  Table  5.1  lists  the  axial 
velocities  calculated  at  the  points  along  the  dendrite  inter¬ 
face.  It  can  be  seen  that  again  all  the  velocities  calculated 
are  within  1%  of  the  theoretical  value. 
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TABLE  5.1 

A  COMPARISON  OF  THE  IVANTSOV'S  METHOD  AND  THE  FINITE  ELEMENT 
METHOD  FOR  THE  AXIAL  GROWTH  VELOCITIES  OF  VARIOUS  POINTS 
ON  THE  INTERFACE  OF  AN  ISOTHERMAL  DENDRITE 
(NON-DIMENSIONAL  TIP  CURVATURE  =  -0.022) 


Position  Along  the 
Interface  (Angle  in 
Degrees  frcm  the 
Z  Axis) 


Growth  Rate  by  Growth  Rate  by 

the  Ivantsov  Method  the  Finite 

(Non-Dimensionalized)  '  Element  Method 

(Non-Dimension- 

alized) 


0 

0.0125 

0.0124 

17 

0.0119 

0.0118 

28 

0.0112 

0.0110 

39 

0.00942 

0.00964 

50 

0.00771 

0.00797 

61 

0.00569 

0.00601 

72 

0.00335 

0.00383 

Axial  Growth  Rate  of  the  Tip 
(Non-Dimensional ) 
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Figure  5.5  A  Comparison  of  the  Finite  Element 

Solution  and  the  Ivantsov  Solution 
for  the  Growth  Rate  of  a  Dendrite  as  a 
Function  of  the  Tip  Radius  of 
Curvature 


' 


Ill 


5.3.2  Non-Isothermal  Dendrite 

Computations  were  next  performed  on  a  non-isothermal 
dendrite  following  the  Bolling  &  Tiller  and  Temkin  models. 
Calculations  were  performed  for  various  values  of  subcoolings 
(i.e.  various  Stefan  numbers  ranging  from  0  to  1).  Although 
the  computer  program  is  based  on  non-dimensionalized  variables, 
the  physical  constants  given  in  Appendix  III  are  used  to 
convert  the  results  to  specific  systems  such  as  ice,  tin, 
succinonitrile  etc.  The  .maximum  velocity  criterion  of  Bolling 
and  Tiller  was  used  to  determine  the  actual  tip  radius  of 
curvature  of  the  dendrite.  The  stability  criterion  of  Langer 
and  Muller-Kriembhaar  was  also  considered  to  obtain  a  compari¬ 
son  with  the  maximum  velocity  principle. 

Figure  5.6  presents  the  non-dimensional  tip  velocity 
of  an  ice  dendrite  as  a  function  of  the  non-dimensional  tip 
curvature  for  various  values  of  the  stefan  number.  It  is 
clear  that  at  each  Stefan  number  the  v~K  curve  passes  through 
a  maximum  which  is  the  actual  tip  curvature  for  a  dendrite 
according  to  the  maximum  velocity  principle.  Figure  5.7 
presents  a  comparison  of  present  calculations  with  those 
obtained  by  the  Temkin  model  and  the  Ivantsov  model  for  a 
Stefan  number  of  0.05.  The  curvature  and  velocity  based  on 
the  L-MK  stability  theory  are  also  shown  in  the  figure. 

Figure  5.8  shows  the  relationship  between  the  actual 
velocity  of  growth  of  ice  dendrites  as  a  function  of  super¬ 
cooling  based  on  available  experimental  correlations  as  well 
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Figure  5.6  Plot  of  the  Growth  Velocity  of  a 

Dendrite-  Growing  in  a  Supercooled 
Water  as  a  Function  of  the  Tip 
Curvature  for  Various  Stefan  Numbers 
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Figure  5.7  Plot  of  Tip  Growth  Velocity  of  a 

Dendrite  as  a  Function  of  the  Tip 

Curvature  k 
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Figure  5.8 


Plot  of  the  Growth  Velocity  of  Ice 
Dendrites  as  a  Function  of  Supercooling 
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as  the  present  calculations  based  on  the  maximum  velocity 
principle.  It  is  seen  that  although  no  two  correlations  agree 
completely ,  the  results  based  on  present  calculations  fall 
well  within  the  range  of  values  obtained  by  different  correla¬ 
tions.  Figure  5.9  compares  the  values  of  tip  radius  of  curva¬ 
ture  obtained  by  the  present  method  with  those  of  previous 
researchers.  Again,  the  values  obtained  here  compare  favour¬ 
ably  with  the  available  data  in  the  literature. 

It  is  to  be  noted  that  all  the  correlations  available 
in  the  literature  show  a  simple  power  law  relationship  between 
the  velocity  of  growth  and  the  amount  of  subcoolings,  i.e. 
all  the  lines  in  Figure  5.8  except  the  one  based  on  the  finite 
element  calculations  have  constant  slope.  The  line  based  on 
the  present  calculations  deviates  from  a  constant  slope  show¬ 
ing  the  true  effect  of  the  varying  Stefan  number. 

Figures  5.10  and  5.11  present  the  comparisons  between 
the  present  method  and  the  available  experimental  data  for 
the  velocity  of  growth  of  tin  and  succinonitrile  dendrites 
respectively.  The  results  show  that  the  present  method  using 
the  maximum  velocity  criterion  gives  growth  velocities  which 
are  consistently  a  factor  of  two  to  three  too  high,  but  are 
consistent  with  other  approximate  analytical  techniques  that 
use  the  maximum  velocity  criterion. 

In  all  of  the  above  calculations  it  was  assumed  that 
the  dendrite  shape  was  a  paraboloid  of  revolution.  Although, 
it  is  agreed  by  most  researchers  that  most  dendrites  are  closer 
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Figure  5.9  Plot  of  the  Tip  Radius  of  Curvature  as 

a  Function  of  Supercooling  for  Ice 
Dendrites 
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Figure  5.10  Plot  of  the  Growth  Velocity  of  Tin 

Dendrites  as  a  Function  of  Supercooling 
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Figure  5.11  Plot  of  the  Growth  Velocity  of 

Succinonitrile  Dendrites  as  a  Function 
of  Supercooling 
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to  a  paraboloid  shape,  many  authors  believe  that  the  actual 
shape  varies  somewhere  between  a  paraboloid  of  revolution  and 
a  parabolic  cylinder.  In  particular  the  two  principle  radii 
of  the  ice  dendrite  appear  to  be  largely  different  [27] . 

Figure  5.12  presents  the  results  based  on  isothermal  and  non- 
isothermal  parabolic  cylinders,  isothermal  and  non-isothermal 
paraboloids,  and  various  experimental  correlations  for  the 
growth  of  ice  dendrites.  It  is  clear  that  the  dendritic  growth 
velocity  will  depend  strongly  on  the  three  dimensional  shape 
of  the  dendrite  as  the  results  based  on  various  models  differ 
by  several  orders  of  magnitude.  Although  the  values  of  velo¬ 
cities  based  on  the  current  non-isothermal  paraboloidal  model 
using  the  maximum  velocity  criterion  come  close  to  the  experi¬ 
mental  data,  this  result  may  be  somewhat  fortuitous. 

Examining  the  results  in  more  detail  will  show  that 
either  the  quasi-steady  state  growth  assumption  model  is  not 
appropriate  for  a  dendrite  or  the  shape  which  obeys  the  quasi¬ 
steady  state  growth  is  other  than  a  paraboloid  of  revolution. 
Holzman  [45] ,  [46]  showed  that  both  the  paraboloid  (needle 

crystal)  and  the  parabolic  cylinder  (platelet)  shapes  do  not 
satisfy  the  heat  balance  conditions.  He  calculated  a  quantity 
he  defined  as  an  excess  velocity  for  both  these  shapes.  He 
showed  that  if  one  calculates  the  velocity  of  interface  along 
each  point,  the  condition  of  steady  growth  is  violated.  Similar 
calculations  by  the  present  method  confirmed  his  results. 

Figures  5.13  and  5.14  show  the  excess  velocity  obtained  by  the 
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Figure  5.12  Plot  of  the  Growth  Velocity  of  Ice 

Dendrites  as  a  Function  of  Supercooling 
for  Various  Dendrite  Tip  Shapes 
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Position  Along  the  Interface 
(Angle  in  Degrees  From  the  Z  Axis) 


Figure  5.13 


Plot  of  the  Growth  Velocity  of  Various 
Points  Along  the  Interface  of  a 
Paraboloidal  Dendrite  k 
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Position  Along  the  Interface 
(Angle  in  Degrees  From  the  Z  Axis) 


Figure  5.14 


Plot  of  the  Growth  Velocity  of  Various 
Points  Along  the  Interface  of  a 
Parabolic  Dendrite  k 
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finite  element  method  for  the  two  shapes  under  consideration. 
It  is  seen  that  the  hypothetical  paraboloid  and  the  parabolic 
cylinder  bulge  behind  the  tip  as  found  by  Holzman.  The  dotted 
line  in  both  figures  represents  the  velocities  required  to 
maintain  the  time  invariant  shape,  while  the  solid  line  shows 
the  actual  velocity  obtained  by  calculation. 

Summarizing  the  results  of  this  chapter,  it  is  eviden¬ 
tly  clear  that  the  finite  element  method  can  be  used  to  study 
the  dendritic  growth  phenomenon.  However,  to  gain  a  better 
understanding  of  dendritic  growth,  it  is  important  to 
simulate  dendrite  growth  over  a  time  period.  The  next  chapter 
is  devoted  to  this  study. 


CHAPTER  6 


SIMULATION  OF  DENDRITIC  GROWTH 

The  purpose  of  this  chapter  is  to  make  an  attempt 
to  study  the  development  of  the  crystal  shape  with  time. 

As  described  previously,  the  experimental  evidence  shows 
that  dendrites  grow  with  a  constant  speed  for  a  given  super¬ 
cooling.  But  the  study  of  dendrite  growth  by  the  finite 
element  model  gave  no  hint  that  the  dendrite  would  grow  with 
a  time  invariant  shape  and  velocity.  The  maximum  velocity 
principle  presented  a  convenient  tool  to  obtain  a  unique 
value  for  the  velocity  of  growth,  but  it  is  not  based  on  any 
strict  physical  reasoning.  The  stability  criterion  of  Langer 
and  Muller-Krumbhaar  has  also  not  been  thoroughly  tested.  The 
aim  here  is  to  make  an  attempt  to  obtain  a  shape  which  is 
time  invariant  and  the  corresponding  velocity  of  growth  of 
a  dendrite. 

To  simulate  dendrite  growth,  the  following  steps  were 
required . 

1.  An  initial  dendrite  shape  was  chosen. 

2.  The  finite  element  program  was  used  to  calculate  heat 

fluxes  at  the  interface  and  from  these  heat  fluxes  the 
velocity  of  advance  at  each  point  on  the  interface  was 
calculated  in  the  same  way  as  it  was  done  in  the  previous 

chapter . 
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3 .  Each  interface  node  was  then  advanced  by  an  amount 
equal  to  the  velocity  times  a  time  step.  This  time  step  was 
chosen  such  that  the  maximum  advance  of  a  node  was  about  0.1 
to  0.2  of  the  tip  radius. 

4.  The  interface  nodes  were  then  shifted  in  the 
negative  Z  direction  so  that  the  finite  element  domain  remained 
centered  on  the  node  at  the  dendrite  tip.  In  this  shift,  the 
nodes  must  be  moved  along  radial  lines  so  that  their  spacing 

is  not  severely  distorted  by  changes  in  dendrite  shape. 

5.  New  interface  curvatures  were  then  calculated  for 
each  nodal  point.  The  technique  used  in  this  step  will  be 
discussed  later. 

6.  The  finite  element  grid  was  then  constructed  again 
using  the  new  interface  position  and  the  velocities  at  the 
interface  nodes  recalculated. 

7.  Steps  (3)  through  (6)  were  then  repeated  a  large 
number  of  times.  The  result  was  a  simulation  of  the  dendrite 
interface  evolution  as  a  function  of  time. 

Perhaps  the  most  difficult  step  in  this  procedure  is 
step  (5) — the  calculation  of  interface  curvature.  Cubic 
splines  were  considered  first  to  represent  the  moving  shape¬ 
changing  interface.  Cubic  splines  were  selected  because  it  is 
known  that  splines  can  fit  various  shapes  easily  and  continuity 
of  slope  and  curvature  is  maintained  at  each  knot.  Ability  to 
calculate  the  curvature  at  each  node  is  essential  since  the 
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curvature  is  needed  to  calculate  the  freezing  point  depression 
due  to  Gibbs— Thompson  effect  at  each  point. 

The  use  of  cubic  splines  was  however,  not  very  success¬ 
ful.  It  was  found  that  the  cubic  spline  routine  interpolated 
the  moving  interface  very  accurately,  but  the  curvatures 
calculated  near  the  tip  were  not  very  satisfactory.  Since 
the  cubic  spline  is  required  to  pass  through  every  point  and 
the  nodal  points  near  the  tip  are  spaced  very  closely,  the 
curvature  values  obtained  fluctuated  wildly  from  one  node  to 
the  next.  This  provided  a  numerical  unstability  that  grew  as 
the  calculation  proceeded. 

To  eliminate  the  above  mentioned  problem  a  slightly 
different  curve-fitting  strategy  was  employed.  Since  the 
difficulty  in  the  curve  fitting  was  near  the  tip,  the  curve 
fitting  routine  was  split  in  two  parts.  Near  the  tip  a  para- 
pola  was  fitted  using  the  least  square  technique  and  for  the 
remaining  part  of  the  interface,  the  interface  curvature  was 
obtained  by  using  the  finite  difference  method  to  evaluate  the 
second  derivative.  In  fitting  a  parabola  to  the  dendrite  tip 
the  node  right  at  the  tip  was  not  included.  As  was  the  case 
for  the  sphere  in  Chapter  4,  the  calculated  velocity  at  this 
point  is  not  consistent  with  that  at  the  adjacent  nodes  pre¬ 
sumably  because  of  the  singular  nature  of  this  point  in  the 
finite  element  grid.  Using  this  technique  a  smoothly  varying 
curvature  that  appeared  reasonable  was  obtained.  However,  the 
calculation  of  curvatures  from  a  set  position  coordinates  for 
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the  interface  remains  as  one  of  the  major  sources  of  potential 
numerical  errors. 

Initially,  to  make  the  long-term  growth  calculations 
simpler,  the  flux  on  the  solid  side  of  the  dendrite  was  neglec¬ 
ted.  In  all  calculations  presented  in  this  section  the  Stefan 
number  was  set  equal  to  zero.  The  results  will,  therefore, 
apply  only  to  cases  where  the  supercooling  is  relatively  small. 
Figures  6.1  through  6.4  present  the  results  obtained  for  the 

long-term  growth  of  dendrites  with  k  =0.  Figure  6.1  shows 

s 

the  evolution  in  shape  of  a  dendrite  interface.  It  is  seen 
that  starting  with  a  paraboloidal  shape,  as  the  dendrite  grows, 
instabilities  develop  near  the  dendrite  tip.  The  tip  becomes 
broader  and  slows  down  before  ultimately  reaching  a  time  invar¬ 
iant  shape  with  constant  velocity  of  growth. 

To  determine  if  the  results  obtained  for  the  long-term 
growth  by  the  finite  element  method  are  dependent  on  the 
initial  conditions,  the  program  was  run  with  two  different 
initial  paraboloids.  Figure  6.2  presents  the  velocity  and 
curvature  obtained  by  using  an  initial  paraboloid  with  broad 
tip  while  Figure  6.3  shows  similar  results  for  an  initial  para¬ 
boloid  with  a  sharp  tip.  The  final  velocity  and  tip  curvature 
for  the  two  initial  shapes  are  within  0.5  percent  of  each  other. 

Figures  6.2  and  6.3  do  show  an  oscillatory  behaviour  as 
would  be  suggested  by  the  Langer— Krumbhaar  stability  analysis; 
however ,  they  also  suggest  that  the  oscillations  are  highly 
damped  so  that  the  dendrite  quickly  approaches  a  time  invariant 
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Figure  6.1  Computer  Plot  of  the  Long  Term  Growth  of  a  Dendrite 

Interface  (k  =  0;  initial  shape  —  Perfect  Paraboloid) 
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Figure  6.2  Plot  of  Growth  Rate  and  Tip  Curvature  Versus  Time  to 

Determine  the  Time  Invariant  Shape  of  a  Dendrite  (k 
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Figure  6.3  Plot  of  Growth  Rate  and  Tip  Curvature  Versus  Time  to 

Determine  the  Time  Invariant  Shape  of  a  Dendrite  (k 
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Figure  6.4  The  Final  Shape  of  a  Time  Invariant  Dendrite  Interface 

and  a  Paroboloid  with  Same  Tip  Curvature 
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shape.  The  ultimate  time-invariant  shape  is  shown  in  Figure 
6.4  and  compared  with  a  paraboloidal  shape  with  the  same  tip 
curvature.  These  results  would  suggest  that  the  dendrite 
obeys  the  quasi-steady  state  assumption,  but  the  shape  is 
slightly  different  than  the  paraboloid  of  revolution  as  assumed 
previously. 

Next,  calculations  were  attempted  including  the  flux 
on  the  solid  side.  With  k  =  4k„  (i.e.  ice-water  system) 

similar  calculations  were  performed  as  described  before. 

Figure  6.5  presents  the  results  of  these  calculations.  Start¬ 
ing  with  two  different  paraboloids,  the  dendrites  reach  a  steady 
velocity  and  tip  curvature. 

The  above  described  results  suggest  that  the  dendrites 
grow  with  a  time-invariant  shape  which  is  slightly  different 
from  a  paraboloid  of  revolution.  Also,  it  is  found  that  the 
velocity  of  growth  of  the  dendrites  is  considerably  lower  than 
that  predicted  by  the  maximum  velocity  principle.  In  fact, 
the  velocity  obtained  by  the  current  method  is  closer  to  that 
predicted  by  the  Langer-Krumbhaar  stability  criterion.  Figure 
6.6  presents  the  comparison  of  velocity  of  growth  for  ice 
dendrites  as  a  function  of  supercooling  obtained  by  different 
methods.  It  is  evident  that  the  velocities  obtained  by  the 
long-term  growth  model  and  the  Langer-Krumbhaar  stability 
criterion  are  very  similar  and  also  in  agreement  with  some 
measurements  for  growth  of  ice  dendrites  made  by  Kallungal  and 
Barduhn  [28]  recently. 


r 


133 


A  xeuoTsuauiTa-uoN 


S200  * 


200  * 


r- 

o 


vO 

o 


m 

o 


o 


m 

o 


CM 

O 


o 

o 


CP 

G 

•H 

■P 

P 

0  0  — 

4->  4-> 

c/> 

0  II 

e  <u 

•h  -p  <n  o* 

H  -r-  44  Lm 

p  w 

W  TJ 
G  G 
in  0  0 
P  Q  P 
0  G 

>  0  -P 

0 

0>w  > 
POP 
G  G 

-P  0  U 

0  0. 

>  0  a 
p  jG  -H 
G  <*)  E-t 
U 

-P  O 

G  *H 

•H  0  1— I 

E-t  -H  O 
P  X) 

0  0  0 
C  >  P 
0  G  0 

W  pL, 

0 

-p  0  -p 

0  e  g 

(X  -H  0 
H  P 
XI  0 

.p  0  iw 

£  x:  <p 

O  -P  *H 
P  Q 

U  0 
C  O 
M-l  -P  S 

o  a  £ 

p 

u  0  E 

O  -P  o 
r-1  0  P 
CU  Q  En 


in 

VO 

0 

P 

G 

CP 

-H 

tu 


AT  °C 


Figure 


6.6 


Plot  of  Dendrite  Growth  Rate  as 
Function  of  Supercooling  AT  for 


a 

ice. 


135 


Although  simulations  were  not  done  for  the  conductivity 

ratio  k  /k^ ,  appropriate  for  tin  or  succinonitrile ,  it  is 

probable  that  the  calculation  of  the  time-invariant  dendrites 

for  these  materials  would  give  velocities  closer  to  what  is 

measured  than  the  use  of  the  maximum  velocity  criterion  would. 

For  k  =0  and  k  /k0  =  4,  the  non-dimensional  velocities  were 
s  s Z 

respectively  a  factor  of  2  and  2.7  less  for  the  time-invariant 
solution  than  the  maximum  velocity  one.  This  is  very  close 
to  the  differences  observed  in  Figures  5.10  and  5.11  between 
the  measured  values  and  the  predictions  based  on  the  maximum 
velocity  assumption. 


CHAPTER  7 


SUMMARY  AND  CONCLUSIONS 

This  work  has  been  an  investigation  of  the  applicabil¬ 
ity  of  the  finite  element  method  to  the  study  of  crystal 
growth  transport  processes.  The  objective  was  to  see  if  the 
finite  method  can  be  applied  to  gain  a  better  understanding 
of  dendritic  growth  phenomenon. 

To  attain  the  objective  set  out  for  this  work,  the 
approach  was  developed  in  three  stages, 

1.  The  first  stage  was  designed  to  demonstrate 

that  the  finite  element  method  can  be  successfully  applied  to 
study  the  problems  of  stability  of  spherical  and  cylindrical 
particles  growing  in  supercooled  liquids.  This  was  done  by 
comparing  the  solution  of  simple  heat  flow  problems  by  the 
analytical  method  with  solution  by  the  finite  element  method. 
Both  approaches  were  shown  to  give  identical  results  (Chapter 
4).  It  was  found  that  the  finite  element  method  could  give 
accurate  results  for  stability  problems  involving  irregular 
boundaries  (perturbed  cylinders  and  spheres). 

2.  In  the  second  stage,  a  finite  element  model  was 
used  to  study  the  dendritic  growth  problem.  This  problem  was 
chosen  because  of  the  great  difficulty  experienced  in  the 
anayltical  treatments  in  the  solution  of  the  diffusion 
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equations  which  govern  the  process  of  dendritic  growth.  The 
choice  was  reinforced  by  the  widespread  occurrence  of  this 
phenomenon ,  which  has  great  technological  importance. 

The  dendritic  growth  problem  was  formulated  such  that 
the  essential  features  of  the  growth  process  were  included 
in  this  model.  A  paraboloid  of  revolution  or  a  parabolic 
cylinder  of  isotropic,  pure  material  was  considered  to  grow 
into  a  pure,  supercooled  melt.  The  heat  flow  in  the  solid 
side  was  also  accounted  for.  The  results  obtained  from  this 
model  were  compared  with  the  available  data  based  on  experi¬ 
mental  as  well  as  analytical  research.  The  results  obtained 
by  the  present  method  agree  very  well  with  the  corresponding 
available  correlations.  However,  the  calculations  from  the 
model  at  this  stage, as  well  as  the  available  analytical  data, 
did  not  fully  explain  the  observed  dendritic  growth  phenomenon. 

3.  In  the  third  stage,  an  attempt  was  made  to  study 
the  long-term  growth  of  a  dendrite.  The  moving  interface  was 
described  in  a  manner  that  could  adequately  represent  the 
dendritic  growth  process.  It  was  determined  that  the  steady 
state  shape  for  a  growing  dendrite  was  close  to  a  paraboloid 
but  not  exactly  a  paraboloid  of  revolution.  It  was  also  found 
that  the  growth  velocity  thus  obtained  was  smaller  than  that 
predicted  by  the  often  used  maximum  velocity  principle.  The 
growth  velocities  were  closer  to  those  predicted  by  the 
Langer-Krumbhaar  stability  criterion. 
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The  most  interesting  future  application  of  the  finite 
element  method  is  to  the  simulation  of  crystal  growth,  such  as 
the  studies  done  in  the  last  chapter.  Although  the  prelimin¬ 
ary  results  presented  in  that  chapter  are  encouraging,  there 
remain  some  major  concerns  about  its  success  in  other  applica¬ 
tions  . 

1.  The  calculation  of  surface  curvatures  are  necessary 
for  the  application  of  the  Gibbs-Thompson  effect  on  interface 
temperatures.  This  can  be  done  with  relative  ease  for  a  sur¬ 
face  described  by  an  analytical  function;  however,  for  a  sur¬ 
face  defined  only  by  a  finite  set  of  coordinate  points  this 
calculation,  which  involves  second  differences  of  the  coordin¬ 
ates,  can  introduce  large  numerical  errors.  Improvements  in  the 
method  by  which  this  calculation  was  done  may  be  required  to 
improve  the  accuracy  with  which  the  method  simulates  the  long¬ 
time  behaviour  of  crystal  growth.  A  method  based  on  coordinate 
transformation  may  reduce  some  of  the  difficulties  faced  here. 

2.  A  second  concern  is  that  numerical  instabilities 
may  be  mistaken  for  or  confused  with  real  instabilities  when 
one  is  using  a  numerical  technique  to  study  a  problem  like 
this.  Generally  one  would  expect  the  numerical  instabilities 
would  be  very  sensitive  to  changes  in  nodal  spacing  whereas  the 
real  instabilities  would  not;  therefore,  repeatability  of  the 
results  for  different  nodal  spacing  is  an  essential  computation¬ 
al  test  of  the  reliability  of  the  predictions. 
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3.  In  simulation  of  crystal  growth  it  is  the  'velocity 
excess  that  determines  the  changes  in  crystal  shape  with  time. 
This  means  that  it  is  not  sufficient  to  calculate  the  mean 
velocity  of,  for  example,  a  dendrite  to  an  acceptable  degree 

of  accuracy,  say  one  percent,  but  one  must  be  concerned  about 
calculating  differences  in  velocities  between  nodal  points  to 
this  accuracy.  In  such  things  as  determining  the  required 
domain  size  or  nodal  spacing  one  should,  therefore,  do  the 
selection  on  the  basis  of  calculated  'velocity  excess'  rather 
than  just  the  mean  growth  rate. 

4.  A  final  concern  is  the  computer  time  required  to 
carry  out  these  simulations.  In  the  calculations  that  were 
shown  in  Chapter  6,  871  nodes  were  normally  used.  This  pro¬ 
duced  21  discrete  coordinate  points  along  the  dendrite  surface 
to  define  its  shape  which  was  considered  to  be  bare  minimum 
needed.  Carrying  out  a  simulation  of  a  dendrite  growth  over 

a  distance  of  the  order  of  30  times  the  tip  radius  required  800 
seconds  of  CPU  time  on  an  Amdahl  460/V7.  One  would  ideally 
like  to  double  the  number  of  nodes  on  -the  dendrite  surface, 
which  would  increase  the  total  number  of  nodes  by  a  factor  of 
two  and  the  computational  time  by  about  a  factor  of  eight. 

Some  runs  were  made  using  this  increased  number  of  nodes;  how¬ 
ever,  the  availability  of  computer  time  became  a  limiting  factor 
on  the  amount  of  work  that  could  be  done.  Obviously  if  extensive 
crystal  growth  simulation  work  was  to  be  done  an  effort  would 
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have  to  be  made  to  reduce  the  computational  time  required  and 
as  well  it  could  only  be  done  on  the  faster  of  the  new  computer 
machines . 

If  these  concerns  and  problems  with  the  use  of  the 
finite  element  method  of  crystal  growth  simulation  can  be  over¬ 
come,  a  number  of  interesting  applications  can  be  envisioned. 

1.  With  a  small  modification  of  the  present  program  the 
time  evolution  of  a  crystal  from  a  spherical  nucleus  could  be 
simulated.  Small  perturbation  analysis  suggests  that  the  Y^0 
mode  is  the  lowest  unstable  mode  of  perturbation  for  a  sphere. 
The  large  amplitude  growth  of  this  perturbation  is,  however, 
unknown. 


2.  Also  with  small  modification  the  effects  of  inter¬ 
face  kinetics  on  the  growth  shapes  can  easily  be  studied.  For 
this  one  just  has  to  introduce  an  additional  term  into  the 
equation  for  the  interface  temperature  that  depends  on  the 
local  interface  velocity. 

3.  With  some  more  extensive  modification  it  should  be 
possible  to  simulate  the  effect  of  solutes  on  crystal  growth. 

This  would  require  the  solution  of  the  mass  diffusion  as  well 
as  heat  diffusion  equations  at  each  time  step. 

4.  Theoretically  it  should  be  possible  to  do  three- 
dimensional  simulations;  however,  at  present  such  work  would 
probably  not  be  feasible  because  of  the  computation  time  required. 
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APPENDIX  1 


THERMODYNAMIC  EQUILIBRIUM  ACROSS  A  CURVED  INTERFACE 

In  the  physical  model  described  in  Chapter  III  it 
was  stated  that  the  equilibrium  fusion  temperature  of  the 
tip  (interface)  of  a  dendrite  is  depressed  an  amount  AT 

r 

owing  to  the  radius  of  curvature  of  the  non-flat  interface. 

The  derivation  for  this  phenomenon,  known  as  capillarity 
(or  Gibbs-Thompson  effect)  will  be  discussed  briefly  in  this 
appendix.  For  a  detailed  study  on  this  subject  the  reader 
is  referred  to  textbooks  on  thermodynamics  such  as  Guggenheim 
[70]  and  Landau  and  Lifshitz  [71]. 

Consider  a  very  small  section  of  the  interface,  over 
which  the  curvature  is  nearly  uniform,  and  with  small  adjacent 
regions  of  the  solid  and  liquid  phases  which  are  nearly 
uniform  in  composition  and  temperature.  A  pressure  difference 
exists  across  the  interface,  which  can  be  determined  by  minimi¬ 
zing  the  free  energy  with  respect  to  an  infinitesimal  change 
of  volume  of  one  phase,  the  total  volume  being  constant.  This 
leads  to  the  condition 


Ps  P Z  YsZ 


d  surface  area 
d  V 


P  and  P  are  the  pressures  at  the  curved  interface,  r  and  r0 
s  Z  i  z 

are  the  principal  radii  of  curvature  and  are  positive  if  the 
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interface  is  concave  toward  the  solid. 


Figure  Al-1  Equilibrium  Across  a  Curved  Interface 


If  the  pressures  P  and  P  are  nearly  the  same  as  the  pressures 

S  X/ 

at  a  flat  interface, 


AP  =  P  -  P  (flat);  |APl  <<  P  . 

s  s  s  s  s 


(A- 2 ) 


it  may  be  appropriate  to  use  the  Gibbs-Duhem  relations, 


Ay 


s 


V  AP 
s  s 


A  T  , 
r 


(A-3) 


S0  A  T  . 
Z  r 


(A- 4 ) 
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At  a  flat  interface  P  (flat)  =  P  (flat) ,  so  from  equation 

^  X/ 


(A-l) 

AP 

s 

-  A  pz  = 

and  Ay  =  Ay  because  the  phases  are  always  in  equilibrium. 

S  X/ 

If  the  liquid  is  assumed  to  be  under  constant  pressure, 
then  AP  =  0,  and  from  equations  (A-2) ,  (A-3) ,  and  (A-4)  it  is 


found  that 

SZ 

-  AT  ( — 
r  \ 

V  y 

-N;>  =  W 

Now,  T f  (S£/N£  - 

S  /N  )  is  the  latent  heat  of  fusion  per 
s  s 

atom,  so 


where  K  =  (—  + 

rl 

Ys£ 

AT  =  -  T-  K,  (A- 7 ) 

r  f  L 

— )  is  the  mean  curvature  of  the  interface. 
r2 
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APPENDIX  2 


DERIVATION  OF  ELEMENT  MATRICES 

In  the  finite  element  formulation  fo  the  dendritic 
growth  phenomenon  described  in  section  3.5.1  the  functional 
minimization  problem  was  reduced  to  a  simple  matrix  equation. 
The  matrix  equation  (3.56)  can  be  simplified  considerably  if 
all  quantities  independent  of  triangle  size  and  shape  are 
evaluated  once  and  for  all.  In  this  appendix  the  derivation 
of  the  element  matrices  for  a  six-node  triangular  element  as 
shown  below  is  briefly  discussed. 


Figure  A2-1  Six  Node  Triangular  Element 
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1.  Matrices  for  Axisymmetric  Formulation 


From  Chapter  III  the  element  matrices  are: 


and 


SA 


m 

ij 


QA 


m 

iJ 


k=l 


k=l 


(P.  PT  +  P.  p_)  dridC 
1  j  l 

( i ,  J  =  1  to  6 ) 


X±Xj  dn  dC 


D 


(A2-1) 


(A2-2) 


with 

<P1,P2,  . .  Pfi>  =-:  (4^-1)  b1/2Am,  (4i|>2-1)  bs/2Am, 

(4if>3-l>  b3/2Am,  2(*2b1  +  iJ;1b2)/Am, 

2(i|/3b2  +  t2t>3)/Am,  2(*1b3  +  tf.3b1)/Am  >  (A2-3) 


<xr  •• 

•  *  •  /  x6>  =  <^1  (2ipx- 

■1)  ,  i|)2  (2i|i2-1)  ,  ^3  ( 2^3-l )  , 

4 *^1^2  ' 

4^2^3 '  4^1 

(A2-4) 

bk  =  "i  -  ^ 

i ,  J,K  cyclic 

(A2-5) 

a„  =  cj  -  a 

i ,  J,K  cyclic 

(A2-6) 
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and 

A”  =  (aib3  -  a3b1)/2 


( A2-7 ) 


The  array  is  found  by  replacing  the  b's  with  a's  in  the 

expression  for  P.. 

1 

To  assist  in  simplifying  the  derivation  of  the  element 
matrices,  define 


*  /* 

*K  (PiPJ  +  PiPJ>  d^C 
J  . 


(A2-8) 


and 


iJ 


=  2 


XixJ  dnd5 


(A2-9) 


Then  in  terms  of  these  quantities 


k=l 


and 


k=l 


(A2-10) 


(A2-11) 


An  explicit  evaluation  of  equations  (A2-8)  and  (A2-9) 
is  obtained  by  a  straight  forward  substitution  of  the  approp¬ 
riate  quantities  into  the  equations  and  making  use  of  Table 

I 

A.  1 .  For  example,  to  evaluate  T^,  the  following  expression 


is  to  be  evaluated: 
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^1  X1X1  dr|d^  (A2-12) 

On  substituting  appropriate  value  of  x-j_  in  terms  of 
the  equation  (A2-12)  becomes 


11 


=  2 


^  (2^-n 2  dndC 


(A2-13 ) 


and  using  the  Table  A.l  for  integration,  results  in 


m 


11 


=  60 


1260 


( A2-14 ) 


Similarly,  all  the  elements  of  the  matrices  R  and  T  can 
be  evaluated  and  thus  the  coefficient  matrices  SA  and  QA 
can  be  obtained.  The  various  matrices  are  listed  below  for 
convenience.  Also,  since  all  these  matrices  are  symmetric, 
only  the  lower  triangular  elements  are  listed.  The  quantity 
a.  is  used  to  represent  the  sum  (a. a  +  b.b  )/60Am. 
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T3 
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60  A' 
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TABLE  A. 1 


COEFFICIENTS  (C)  FOR  AREA  INTEGRALS  IN  AREA- 
COORDINATE  SYSTEM  (FROM  REFERENCE  [  6 ] ) 


Order 

n  =  ei+Bj+6k 

Bi 

6j 

6k 

c 

1 

1 

0 

0 

1/3 

2 

2 

0 

0 

2/12 

1 

1 

0 

1/12 

3 

3 

0 

0 

1/10 

- 

2 

1 

0 

2/60 

1 

1 

1 

1/60 

Remark : 


/*  /• 


and  i,J,k  represent  any  permutation  of 


1/2,3. 


. 
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2-  Matrices  for  Two-Dimensional  Formulation 


are  : 


and 


The  element  matrices  for  two-dimensional  formulation 


.m 

’iJ 


QmT 
1 J 


m 


(P.  P  +  P  P  )  dA  (i , J=1  to  6)  (A2-15 ) 


X .  Xt 
Ai  AJ 


(A2-16 ) 


m 


with  arrays  Pi ,  P  X  Xj  and  Am  defined  exactly  the  same 
as  those  in  previous  section,  except  that  the  coordinates 

n /  £,  now  represent  two-dimensional  system.  The  integrals 
can  be  evaluated  by  following  the  same  procedures  described 
previously  for  axisymmetric  case,  and  one  can  obtain  the 
corresponding  matrices  for  two-dimensional  formulation,  which 
are  listed  below.  Here  or  j 


m 


(a . a_  +  b . b  ) /12A  . 
l  J  l  j 


is  used  to  represent  the  sum 
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and 


QiJ  1260  6 

1  6 
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APPENDIX  3 


AUTOMATIC  MESH  GENERATION 

As  mentioned  in  chapter  3 ,  the  element  mesh  is 
generated  automatically  by  the  computer  on  supplying  a  minimum 
amount  of  data.  For  a  domain  which  is  divided  into  twenty 
blocks  similar  to  one  shown  in  Figure  3.9,  the  following 
input  data  is  required: 

1.  The  z  coordinate  of  nodes  1  and  21. 

2 .  The  angle  between  the  line  drawn  from  the  origin 

to  the  node  and  the  z  axis  for  nodes  43,  63,  85,  105,  127,  147, 
. ,  799,  819. 

3.  The  desired  spacing  between  the  nodes  1,  3,  5,  7, 

. .  21  as  fraction  of  the  total  distance  between  the  nodes 

1  and  21. 

With  the  above  information,  the  r  and  z  coordinates 

of  nodes  43,  63,  85,  105,  127,  147,  . .  799,  819,  841, 

and  861  can  be  obtained  provided  the  shape  of  the  interface 
and  the  outer  boundary  is  known  (i.e.  spherical,  parabolic 
etc.).  The  coordinates  for  remaining  nodes  can  be  determined 
by  using  the  fractional  ratios  of  distances  between  the  nodes 

1  and  21,  4  3  and  63  ,  8  5  and  125,  . .  841  and  8  61.  The  node 

numbering  on  the  solid  side  is  similar  to  the  numbering  on  the 
liquid  side.  The  interface  node  members  are  common  between 
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the  solid  and  liquid  sides.  Such  a  mesh  generation  scheme 
facilitates  quick  and  easy  change  in  spacing  of  nodes  as 
required  by  the  problem  under  investigation. 
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APPENDIX  4 


PHYSICAL  CONSTANTS  FOR  WATER ,  TIN ,  AND  SUCCINONITRILE 

USED  IN  THIS  WORK 


( 1 )  WATER 

Equilibrium  melting  temperature 
Latent  heat  of  fusion 
Density  of  solid 
Density  of  liquid 
Thermal  conductivity  of  liquid 
Thermal  conductivity  of  solid 
Specific  heat  of  liquid 
Thermal  diffusivity  of  liquid 
Interfacial  surface  energy 


273  K 

3 

80  cal/cm 

3 

0.91  gm/cm 
1  gm/cm^ 

0.00133  cal/cm. deg . sec 

0.00532  cal/cm . deg . sec 

1  cal/gm  °C 
2 

0.00144  cm  /sec 
2 

7.17  cal/cm 


(2)  TIN 

Equilibrium  melting  temperature 
Latent  heat  of  fusion 
Density  of  solid 
Density  of  liquid 
Thermal  conductivity  of  solid 
Thermal  conductivity  of  liquid 
Specific  heat  of  liquid 
Thermal  diffusivity  of  liquid 
Interfacial  surface  energy 


505  K 

102  cal/cm^ 

3 

7 . 3  gm/cm 

6.98  gm/cm^ 

0.144  cal/gm. deg . sec 

0.064  cal/cm. deg . sec 

0.060  cal/gm  °C 
2 

0 . 2  cm  /sec . 
13.03  x  10  ^  cal/cm^ 
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SUCCINONITRILE 

Equilibrium  melting  temperature 

331.24  K 

Latent  heat  of  fusion 

3 

10.72  cal/cm 

Density  of  solid 

3 

1.016  gm/cm 

Density  of  liquid 

3 

0.970  gm/cm 

Thermal  conductivity  of  solid 

5.36 

-4 

x  10  cal/cm . deg . sec 

Thermal  conductivity  of  liquid 

5.32 

-4 

x  10  cal/cm. deg . sec 

Specific  heat  of  liquid 

0.478  cal/gm  °C 

Thermal  diffusivity  of  liquid 

2 

0.00116  cm  /sec . 

Interfacial  surface  energy 

2 . 14xl0~7  cal/cm2 

• 

